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ABSTRACT 


A second-order epsilon method is developed for trajec- 
tory optimization problems. The method is applied to sever- 
al aircraft and missile performance and air combat maneuver- 
ing problems. Heavy emphasis is placed on the realistic 
modeling of the flight vehicle's motion and maneuvering 
timitations. 

The proposed optimization technique, which is an exten- 
sion of Balakrishnan's epsilon method, uses either the full 
second-order Newton-Raphson method or the "modified" Newton- 
Raphson method to Minimize the epsilon functional. The full 
Newton-Raphson method exhibits terminal convergence charac- 
teristics superior to the "modified" method, whereas the 
"modified" method is generally superior in the initial 
stages of a problem. An algorithm is developed which uses 
both techniques in a complementary way. 

A new penalty functional which has desirable theoretical 
properties and exhibits excellent computational behavior is 


introduced to treat state and control inequality constraints. 
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I. INTRODUCTION 


The objective of the research reported on herein is to 
develop a method of solving realistic problems in aircraft 
and missile performance optimization. Optimization problems 
of this type have been the subject of considerable research 
Reis, 45. 2, 3, 4057-5546, and 7]. “he mathematical models 
used in these references are the products of many simplifi- 
cations and assumptions. Typically, the degree of simplifi- 
cation used to render these problems solvable by some opti- 
mization technique is such that the solutions obtained are 
of limited practical value. This is particularly true in 
the modeling of aircraft maneuvering limitations, such as 
aerodynamic stall, maximum structural load factor, and 
Macará Mach number, which require the use of multiple 
state and control inequality constraints. Since these 
limitations play an extremely important role in maneuvering 
flight and air combat, they must be modeled as accurately 
Sie ible. Tt is, therefore, imperative that the optimi- 
Zation technique used to solve the problems posed herein 
be capable of handling state and control inequality 
foodies With relative ease. 

Balakrishnan's epsilon method is an attractive optimiza- 
tion technique because of the natural manner in which state 
and control inequality constraints are introduced. The 


epsilon method is a penalty method in which terms are added 








to the performance measure to penalize deviations from the 
state equations written as equality constraints. Likewise, 
state and control inequality constraints may be treated by 
the addition of appropriate penalty terms to the performance 
measure. The resulting augmented performance measure is 
minimized by an appropriate algorithm for solving unconstrained 
optimization problens. I 
The optimization technique used most successfully in 
the literature with the epsilon method is a "modified!" 
Newton-Raphson technique, hereafter referred to as the MNR 
technique. In this method certain second-order terms pre- 
sent in the full Newton-Raphson formulation are neglected. 
It is argued [Ref. 8] that since computer storage and time 
requirements to compute these terms are large, and since 
satisfactory results can be obtained over a large class of 
problems without the terms, their inclusion is not justified. 
For these reasons the full Newton-Raphson formulation, here- 
after referred to as the FNR technique, has not been 
previously utilized with the epsilon method. 
Difficulties were experienced by the author, however, 
in applying the MNR technique to problems of the type formu- 
lated in this dissertation. The MNR technique was not 
effective in problems with state equations and multiple 
inequality constraints resulting from a realistic modeling 
of the flight vehicle's motion and maneuvering limitations. 
For this reason the FNR method was investigated and 


found to be feasible in terms of computational storage and 
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time requirements. The FNR method exhibits terminal con- 
vergence characteristics superior to the MNR method although 
the MNR method is generally superior in starting a problem. 
Problems not solvable by the MNR method alone were solved 
by an algorithm which uses both methods in a complementary 
way. The power of the FNR technique close to the minimum 
can also be used to advantage to obtain a family of optimal 
trajectories for different end conditions. The optimal 
trajectory for one set of end conditions is used as a first 
guess for the optimal trajectory for a neighboring set of 
end conditions. 

Several simplified problems in aircraft performance 
optimization were attempted initially to gain experience 
with the epsilon method. In these problems inequality 
constraints were treated by using interior penalty func- 
tionals of the type recommended by Fiacco, McCormick, and 
Jones [Refs. 9 and 10]. Computational results were unsatis- 
factory. Difficulties were experienced in keeping the con- 
strained state or control completely admissible; a require- 
ment for the e of an algorithm with this type of 
penalty functional. To alleviate this difficulty a new 
penalty functional for inequality constraints is introduced 
which exhibits excellent computational behavior. The pro- 
posed functional has performed well in computation with up 
to eight inequality constraints represented in a single 


Problem. 





The thesis is divided into eight sections. In Section 
II the epsilon method is presented. The FNR and MNR tech- 
niques are derived and discussed. The effectiveness of 
the FNR method as opposed to the MNR method is demonstrated 
by a scalar example. Finally, the computational experience 
galned with both methods in solving realistic performance 
problems is presented. In Section III the method of treating 
state and control inequality constraints is presented. The 
author's experience with interior penalty methods is related 
and a new penalty functional is proposed. Computational 
experience with the new penalty functional is related and, 
finally, several desirable theoretical properties of the 
new penalty functional are presented. In Section IV the 
algorithm developed for minimizing the epsilon functional 
by either the MNR or FNR methods is presented. In Sections 
V, VI, and VII three aircraft and missile performance opti- 
mization problems are solved. These problems are pertinent 
and realistic in their operational applicability. The 
three-degree-of-freedom models are the same as those used 
in basic aircraft performance analysis. Finally, the 


summary and conclusions are presented in Section VIII. 





II. THE EPSILON METHOD 


This section describes the epsilon method and reviews 
the significant contributions of other investigators. The 
full Newton-Raphson (FNR) equations for minimizing the aug- 
mented performance measure are derived and compared to the 
modified Newton-Raphson (MNR) equations published elsewhere 


[Refs. 8, 11, and 12]. 


A. DESCRIPTION OF THE EPSILON METHOD 
l. Statement of the Problem 
A dynamic system characterized by the nonlinear 


state equations 
x(t) = f[x(t),u(t),t] (2.1) 
is to be controlled to minimize the performance measure 
T 
J(x,u) = h[x(T),T] + f gLx(t) ,u(t),tJat (2.2) 
an = t x = 
O 


where x(t) is an n x 1 state vector and u(t) is an £ x 1 
control vector. State and control inequality constraints 
are omitted for the present. In Section III the inclusion 
of these constraints is discussed in detail. 
2. The Augmented Performance Measure 
In the epsilon method as proposed by Balakrishnan 


(Refs. 13 and 14], the performance measure (2.2) is augmented 
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by a penalty functional which involves a weighted integral 
of the Euclidean norm of the state equations written as 


equality constraints. The augmented performance measure is 


l 
|= 


T 2 
or) = JO um f ||x - f(x,u,t)|| at (o 
E 2 E > 


J(x,u) + 


m | 


3. (x,u). (2.4) 


The weighting ee e is a positive quantity. 

3. Behavior as € » 0 

As e is reduced, the penalty term d is more heavily 

weighted, thereby placing greater emphasis on satisfying 
the state equations. Balakrishnan [Refs. 13 and 14] and 
Taylor [Ref. 11] have shown that under appropriate assump- 
tions as € + 0, the epsilon method yields the necessary 
conditions of optimality obtained by applying Pontryagin's 


minimum principle [Refs. 15 and 16]: 


k*(6) = gp [Z*E už (t),t], a 
pe(t) = = Sk DX8(t).us(e),tJ, (2.6) 

and 
HLx*(t) u*(t),p*(t),t) « HDx*(t),u(t),p*(t),t] — (2.7) 


where H is the Hamiltonian function defined as 


He 


> A 








HEx(t),u(t),p(t),t] = glx(t),u(t),t] + pi (t)flx(t) ult) t], 


(2.8) 


p(t) is the costate or adjoint vector, u#(t) is an extremal 
control vector, and x*(t) is an extremal trajectory. The 
assumptions made are that the minimization problem has a 
unique solution with x(t) absolutely continuous for each €, 
and that f and g are continuously differentiable in X and 
u [Ref. 11]. Thus, under appropriate assumptions, it can 
be shown that as g > 0, the epsilon method yields the 
results of Pontryagin's minimum principle. That is, if the 
optimal control u*(t,e) of equation (2.3) exists for each 
€, that solution will approach the optimal control u*(t) 
of equation (2.2) as £ > 0. | 

4. State Equation Integration 

lt should be noted that the epsilon method is a 

non-dynamic method in that the state equations are not 
integrated during the minimization process. Once the aug- 
mented performance measure has been minimized a check on 
the degree of satisfaction of the state equations can be 
obtained by integrating the state equations with the optimal 


control. 


` 


B. MINIMIZING THE AUGMENTED PERFORMANCE MEASURE 


i. Sequence of Unconstrained Problems 


Once the augmented performance measure (2.3) is 


formulated, any unconstrained optimization algorithm can 
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be applied to it. A sequence of unconstrained problens 


referred to as 


sub-problems is solved. In each sub-problem 


e is held constant and a minimization is performed until 


some stopping criterion is satisfied. At this point e is 


decreased and a new sub-problem is commenced using the 


optimum trajectory found in the previous sub-problem as a 


first guess. 


In this manner a sequence of sub-problems is 


solved until, 1f convergence occurs, some overall stopping 


criterion is satisfied. 


2. Unknowns and Time Discretization 


The states and controls can be approximated by any 


orthogonal expansions. The coefficients in these expansions, 


along with all 
or unknowns in 
the form (2.9) 
the period can 
Siero at” tne 


time-invariant 


free end conditions, become the parameters 

the optimization. A functional expansion of 
is convenient because it is continuous and 

be selected so that the value of the expansion 
end points. Since the problems solved involve 


systems, E, is selected as zero and the states 


and controls are written as 


x(t) 
age) =| === 
u(t) 


where M is the 


sin — 


= (0) + y(T) - y(0) 4 4 p | sin = (2.9) 


t 


1 T 


Sin Mit 


number of harmonics used and 
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~X 
psa (2.10) 


~U 


u 


is an (n + g) x M matrix of coefficients. The derivative 
of the expansion of the state vector given in equation (2.9) 


with respect to time is required ..and is given by 


COS T 


x(T) - x(0) 


x(t) = = + D 


~X 


T 
T 
am ent 
T x|T 


COS IT (2.11) 


° 


MT Mrt 
T COS T 


The objective is to find the D matrix along with the values 
of the free end conditions which minimize equation (2.3) 
for a given e. In order to perform this minimization, the 
time interval T is divided into (K - 1) sub-intervals each 
of duration At so that there are K discrete time points. 
The augmented performance measure given by Equation (2.3) 


is written as 


II 


T 2 T 2 
1 ; , 
J (x,u,e) uf [xi Ct)-f4 (x,u)] dt + f, E re at 


(2712 


+ 


T 2 
... + f, Baar) | at 


T 
+ J elx,ulät * h[x(T),T]. 
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Suppressing the arguments for clarity, equation (2.12) is 
expanded to yield : 
St 
E gc 
elf (xf) dt 


2At | > KA > 
+f (x,-f) dt + °° ° +f (xi -f4) dt 
At (K-1)At 


= 


T M Mn c E 
f (x,-f,) dt +f (x,-f,) dt + ee. +f (x,-f,) at 
0 At (K-1)4t 


e (2.1) 


At > At > KAt ` 5 
+ 7 (x 7f.) dt + (x -f) dt + --- + de O ae 
YO. At (K-1) At 


t At KAt 
+ T gt Ha ge se f gne nDCD,T] 
0 (K-1)At a 


At 


which can be approximated by 


: 2 
Ja 3 (0) - £, Ex(0) ,u(0) J} " + 4x, (at) - £, Ex(at) u(t) I} - 
Bo {4 (CK-L at) - f [xCIK-1At) su (K-14) 7) a 


° 2 E ’ 


I : aU 
eet caeno - t ponas au e2201) Se ie t (2.14) 


2 
Ix (0) m f E00) ,(0))) 4 {ix (at) = £ (x(t) u(at)} = 
Boe fx. (K-1)4t) - f, E O28) „ut (Kar) M 


fel x (K-DAt) netas + alamo. 
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Hence, the augmented performance measure can be written as 
J. = w.“ + w.“ +... + Wo 2.15) 
= wow (2.16) 


where w is a Q x 1 column vector. The first K elements of 


ware 


w = (* [(=1)6t] z £0x((k-1) at) al nane] | 


(2210) 
Ieee ae 


etc. The form of equation (2.16) is convenient for computer 
programming the epsilon method and for the derivation of 
the minimization techniques that follow. In minimum time 


problems where the performance measure is given by 


T 
J -f dt. (2.18) 
0 


one element of W of the form 


1 
wars [e u sc]: (2.19) 
is used to represent equation (2.18). In these type problems 
the number of time points K is held constant during the 
minimization in order to keep the dimensions of all vectors 


and matrices constant and the time interval At is minimized. 
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The values of the states and controls required in 
equation (2.14) are obtained by evaluating equations (2.9) 
and (2.11) at each time point t = (k - 1)At where k = 1,2, 


...¿K. Written in discrete form equation (2.9) is 


35 T(k-1) 


E E 


y L(K-1)4t]-y(0) "m 
y[(k-1)At] = y(0) + £— — ———— (k-1) * D K-1 


K-1 
š Mr(k-1) 
sin —— 
(2.20) 
and equation (2.11) is 
T de T(k-1) 
K-1)At K-1 
: x[ (K-1)4t ]-x(0) 2T cos 2T(k-1) 
x[(k-1)At] = uc -— a= Da (K-1) At K-1 


(K-1) At 


Mr Mr(k-1) 
Beat ra 


01) 


A vector of unknowns c is formed and is given by 


T - ° ° © 
AA CAG ass 
date, 1? dni, 2? inae dA, M? Zi Z5» (45 Zp> At) 
(2.22) 


i i 





where di 5 is the element in the g bh row and gen column 
of D and Z1» Zoo see Zp represent B oree cnd conditions 
some of which occur at t = 0, and others at t = T. Some 
of the z's correspond to states and others to controls. 


The last element, At, is present only if time is to be 


minimized. The ¢ vector consists of L elements where 

L = (n+%&) x M + P (2.23) 
for all problems except minimum time problems and 

L = (n+g) x M + P + 1 (2.21) 


for minimum-time problems. 

With the states and controls given by equation 
(2.20) and the augmented performance measure given by 
equation (2.16), the problem has been transformed into a 
parameter optimization problem with the unknowns given Dy 
c (2.22). 

3E Minimization Techniques 

The methods which have received attention in the 
literature for finding es which minimizes the augmented 
performance measure given in equation (2.3) are the gradient 
method and a "modified" Newton-Raphson method (MNR). 

The gradient method has been investigated by J. 
Taylor [Refs. 11 and 12] and L. Taylor [Ref. 8] with 


unsatisfactory results. These investigators report that 
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in non-linear problems the gradient method frequently 
obtains false minima and requires considerable computation 
time compared to other methods. 

An MNR method in which certain second-order terms 
present in the full Newton-Raphson (FNR) method are 
neglected has enjoyed greater success and requires less 
computation time than the gradient method [Refs. 8, 11, and 
12]. However, in Ref. 8 difficulties are reported with the 
MNR method in non-linear problems. J, often begins an 
oscillation after two or three iterations and does not 
settle to a minimum. In the problems solved herein the 
same oscillations have been observed when the MNR algorithm 
has been used. Convergence to the minimum, when it does 
occur, is typically very slow. Typical performance of the 


MNR method is shown in Figure 1. 


divergence 


oscillation 


slow convergence 


Augmented Performance Measure 





Iteration number 


Figure Ten 


Augmented performance measure vs. iteration 
number-MNR method 
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The FNR method has not been used by other investi- 
gators with the epsilon method because of the increased 
computer time for each iteration, the additional storage 
space required, and a significantly increased analytic 
workload involved in deriving second partial derivatives. 
Because of the poor performance of the MNR method on problems 
of the type solved herein, the FNR method has been investi- 
gated in detail in this work. With careful programming the 
computer time for each iteration and storage space required 
Tor the FNR method has been reduced to an extent which makes 
the method computationally feasible. 

H. The Full and Modified Newton-Raphson Equations 

The FNR equations for finding c* are derived here 
in a manner which permits the MNR equations derived in the 
literature [Refs. 8 and 11] to be obtained by neglecting 
DOT. 

The augmented performance measure is expanded in a 
Taylor series including up to second-order terms and is 


written as 


t Y i T 2 
J (ctAc) x Is (V oe Ac + %(Ac) iue ey oS (2.25) 


where 


20 | 














_ a 
907 
od, 
90, 
mr Ë 
QJ 
_ A 
9c, 
and 
2 2 2 
9 J. 9 Ja š 9 
2 | 
9°, 9 C4 9C, °C 
97 923 3* 
a a 
2 
96-90 ac IC-9C 
g ey E 2 `l 2 
c a 
34J 34J 3^ 
a a 
9C, 904 9C, 9C. IC 


Oving for the increment of Ja we have 


A 
AJ _ = J,(c+Ac) - J,(c) 


= T I T,„ 2 
SUV eds Ac + (Ac) e J,), (4e). 


C 


al 


(2.26) 


(2.20) 


(2:20) 





Applying the necessary condition for a minimum, we have 


W) , 
Mic) |. Gele j va IR eS (2230 


which when solved for Ac yields 
Ac = A(wés y)? (3) (2.31) 
3 c ace a es i 


If the augmented performance measure is written as 


J Š wÍ (2.14) 


Ux 
I= 


where w is a Q-dimensional column vector, then 


T 
VET = V (ww) 
we Wr: (n 
= 2 (V w) Tw 
and 
2 E" T 
Ve ee | 
~ AUS (2.33) 


2 [Gu (un + (v w) w]. 


as ns 


The matrix N is given by 


ns 
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OW, 9 W. 9W4 
9 C4 9 C, IC, 
dW. 9 W. 9 W. 
A oc ac ac 
Vow = i E : (2.31) 
9 WQ d Wo Wo 
aC, 9C, ð Cr 


2 . a ° ° 
ME w is a three-dimensional array composed of L matrices each 




















š 32w 
of dimension Q x L which has as its ijk th element I 5 
I 
that is 
2 2 2 
0 Wa 0 Wy 0 Wy 
2 
9 C4 90,90, 904907, 
2 2 2 
9 Wo 3 Wo 9 Wo 9 Wa 
2 
9C4 9C, 9C, dC, dC) oc 9c, 
2 
p 9 Wo 
9C,09Cr 
2 a Wo awe a Wo A 
pw - 2 — t NO. — — ie : ° 2 (2.35) 
2 9C4 9C4 9C, ðC} Cr 127 
2 2 2 2 
. 19 
0 Wo 0 Wo e 0 Wo . Wo 
2 h 2 
90,904 IC, ICC 9er 
2 2 DN 
o W ð W ò W 
Q AS Lu. E 
dC, dC, 9C, 9C, 9c, 
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Substituting equations (2.32) and (2.33) into equation (2.31), 


we have 


C 


f.p Pup Ap f.p ^ 


be =-[(VQW) on, * (uou egenus s (36 


This is the full Newton-Raphson equation. The modified 

Newton-Raphson equation can be obtained by neglecting the 

second term in the inverse in equation (2.36), which yields 
ET. T -1 T 

ag [NUI | (a) 7 w (2.37) 


~ ^C 


~~ f^ ^ Ap Ap e ^ 


Several Een concerning equations (2.36) and (2.37) are 
Bm order: 
a. the term VW given by equation (2.34) is a Q x L matrix; 
b. the term (Vw) (Vow) is a symmetric L x L matrix; 
c. the term (V w) w E an L x 1 vector; 
d. the term Y ¿EW given by equation (2.35) is a Q x L x L 
CN rsionamn arpay; 
e. the term (v, ^w) ^w is a symmetric L x L matrix; 


E is defined as the three-dimensional array 


f. (v ^w) 
obtained by transposing each individual Q x L matrix 
given in equation (2.35); 

g. the result of the operation (V ^w) Tw is defined as 


th 


an L x L matrix in which the i column is the 


9 it 
product [5 Men) W. 
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5. A Scalar Illustration of the MNR and FNR Methods 
The potential importance of the second-order term 
neglected in the MNR equations can be illustrated with a 
simple scalar problem in function minimization. The Newton- 
Raphson equation to minimize a function f(x) takes the well 


known form 


p 


Ax = - fix) : (2. 38) 
T£ 
f(x) = w°(x) (2.39) 


which is the form of equation (2.16) with w taken as a 


scalar, then 


s; e (x) (2.40) 


PUR) 


and 


f" ( x) 


2 
2 ES Go | ES — (COR (2.41) 


The FNR equation (2.36) is 


-w(x) $* (x) 
(2.42) 





Ax = 


e 
dw e d w 
| &x j] * w60 a) 


whereas the MNR equation (2.37) is 
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dx 
Ax = a (2.43) 
W 
w(x) 
dw (x) (2.44) 
dx 
Applying equations (2.42) and (2.44) to the function 
MO E eT (2.45) 
in which 
w(x) = (+ m (2.46) 
we obtain 
Ax = - = | 
X = = E (x-1) qo 
for the FNR algorithm and 
CE + 1 
Nenn ST a (2.48) 


2(x-1)? 


for the MNR algorithm. Tables 1 and 2 show the first few 


iterations by both methods from an initial guess of x(0) = 3. 
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MNR Equation (2.48) 





1 3.000 l 17.000 
2 1.937 ee 
3 0.862 . 1.000 
T 190.932 1.329 x 10? 


Table 1 


FNR Equation (2.47) 


Reece [= [= 





: 

3 er 
H c 
p die 
6 12 
T 1% 
8 T: 
9 ! le 


Table 2 


Clearly the MNR equation causes x to diverge after an initial 
period of convergence while the FNR equation causes x to 


approach the minimum. 


6. Computation Experience With the FNR and MNR Methods 
The performance of the MNR method in the preceding 


scalar example is typical of the performance observed by the 


author in large problems. However, the FNR equation is also 
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not uniformly effective when used exclusively in large 
problems. Fortunately, the areas of effectiveness of the 
two methods are complementary. 

In order to discuss the effectiveness of the two 
methods it is convenient to define two areas in the minimiza- 
tion process. Initial behavior refers to the behavior of J, 
during the first two or three iterations in a sub-problem. 
Terminal behavior refers to the behavior of J, after the 
first two or three iterations within the same sub-problen. 
The following behavior has been observed. 

a. Initial behavior: The MNR equation outperforms 
the FNR equation in this area. The ability of the FNR 
equation to minimize Ja is very sensitive to the starting 
value of the unknowns (c). With the values of c far 
removed from the optimum, the FNR equation generally causes 
J, to increase rapidly and diverge from the minimum. ‘The 
MNR equation on the other hand is relatively insensitive to 
the starting c and ean usually be counted on to move J, 
toward the minimum for at least one or two iterations. 

b. Terminal behavior: As the minimum is approached, 
the MNR equation produces the behavior shown in Figure 1. 

The FNR equation, however, generally becomes extremely 
effective in rapidly finding the minimum. 
T. A Combination FNR-MNR Minimization Method 

The obvious approach suggested by the previous obser- 

vations is to devise an algorithm which minimizes by the MNR 


equation initially in a given sub-problem and switches to 
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the FNR equation at some appropriate point in the iteration 


process. Such an algorithm is presented in Section IV. 
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III. AN INEQUALITY CONSTRAINT PENALTY FUNCTIONAL 


In this section a new penalty functional is introduced 


for state and control inequality constraints of the form 


ps 
A 
A 
ps 
` 
pe 
lI 
= 
` 
DD 
` 
w 
H 


4 Í x (t) < 


|^ 


n , te[t ST er, 


c — > 


ND jy 


f 
^ 


< u 


A 
= 
` 
C 
lI 
ER 
` 
no 
` 
` 
HH 
A 


te[t.,T] . (3.2) 


All state and control inequality constraints entered in 
the problems solved herein are of this type. The difficul- 
ties encountered with existing penalty methods which led to 
the use of a new functional are related. The new penalty 
functional has performed well in computation and is used 
exclusively in the solution of the problems presented. 
Additionally, several desirable theoretical properties of 


the proposed penalty functional are presented. 


A. INTERIOR PENALTY METHODS 
1. Past Research 
In Ref. 10 Jones and McCormick present a number of 
theoretical results concerning interior penalty functionals 
of the Fiacco-McCormick type [Ref. 9] in conjunction with 
the epsilon method. If, for example, a state or control, 


denoted by y(t) for generality, is constrained by 


y(t) < Y , telt T] > (3.3) 
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a Fiacco-McCormick penalty functional [Ref. 9] of the form 


Eo 


T 
f —— st | (3.4) 


is added to the augmented performance measure. The behavior 
of the integrand of expression (3.4) for a fixed time 
t e [t TJ as the positive weighting factor r approaches 0 


is shown in Figure 2. 





decreasing r 


L 
Fe 


Figure 2 


Fiacco-McCormick penalty function vs. constrained variable 
for a fixed time 
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If penalty functionals of this type are added to the perform- 
ance measure, it can be shown [Ref. 10] that as r approaches 
0 and € approaches 0, the epsilon method yields Pontryagin's 
minimum principle. The development parallels and augments 
Balakrishnan's work [Refs. 13 and 14] without inequality 
constraints. No computational results, however, are 
presented. 
2. Computational Experience 

A simple problem involving one state variable and 
one constrained control was attempted using the epsilon 
method and a penalty functional of the mes given by 
equation (3.4). The optimal control was on the constraint 
boundary. The algorithm was unable to solve the problem 
by either the FNR or MNR method from a variety of starting 
points. Once the control penetrated the constraint boundary 
for a finite time interval, the algorithm failed on the next 
iteration. The value of r required to keep the control 
admissible for all telt T] throughout the iteration 
process was large, resulting in the augmented performance 
measure being dominated by the Fiacco-McCormick penalty 
term. As a result, the optimal solution [u*(t,e,r)] to the 
augmented problem could not be made to approach the optimal 


solution Lu*(t)]. 


B. A NEW PENALTY FUNCTIONAL: COMPUTATIONAL PROPERTIES 


1. The Form of the New Penalty Functional 


Consider a control or state y(t) which is subject to 


a constraint of the form 
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y(t) e Lys Y y? : telt T] . (3.5) 


A penalty functional of the form 


Tf2a2y(t) -— Ym - y, 72K 
TN) sss f u P at (3.6) 
O 
where Ko is a positive integer ís added to the augmented 
performance measure. The effect of K. can be seen from 


Figure 3 which shows the integrand of equation (3.6) as a 


function of y(t) for a fixed time t e[t >T]. 


increasing Ko 


pu Y 


NA 


Figure 3 


New penalty function vs. constrained variable for a fixed time 
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A functional of the form given by equation (3.6) is added 

to the augmented performance measure for each inequality 
constraint of the form given by equations (3.1) and (3,2). 
For I. control constraints and I, state constraints of this 
form the total augmented performance measure for the epsilon 


method written for time invariant problems with t. = 0 is 


7 2 
nti = [Ix - f(x, ]| dt 


(3.7) 
T E ae Srey. ss 
lee. Ser). 
=) "dis Jr, i71 m KIM 
1 
EU. t1 (3.8) 


The "two-sided" feature of the penalty functional makes it 
especially suited to constraints of the form given by 
equations (3.1) and (3.2). In effect, two inequality 
constraints are included in one penalty term. 
2. Computational Strategy 

The power is increased —Á in numerical 
computation in the same manner as e is decreased. Thus, 
increasingly refined boundaries to the admissible region 
are provided. Both e and > are held constant within a 
sub-problem and are altered between sub-problems. The 
weighting factor r, which is required to provide an overall 
weighting among J, Jas and Js is held constant throughout 


the entire problem. 
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Computational results with this penalty functional 
have been excellent. Up to eight inequality constraints have 


been treated successfully in one problem. 


C. THEORETICAL PROPERTIES OF THE NEW PENALTY FUNCTIONAL 
m = Tntrodtetton 
Three desirable properties of the new penalty 
functional are presented here. These properties and their 
importance are discussed followed by a proof of each 
property. 
a. Penalty functionals of the form of equation (3.6) 


are convex on RP where Cy E R? is defined by 


T 


II > 


y [ay a... .. ay y(t ),y(T)] (3.9) 
and 
y(D-y(t.) mn(t-t_) 
y(t) = y(t.) + = (t-t ) + >. sin m (3.10) 
O m~ O 


It is desirable that the augmented performance measure be 
convex in the unknowns of the minimization to insure that a 
global minimum is attained. If it can be shown that the 
inequality constraint penalty functionals (3.6) are convex, 
then the addition of any number of these functionals does 
not destroy a convexity condition which exists without these 
terms, because the sum of convex functionals is also convex. 
Indeed, me addition of terms of the type given in equation 
(3.6) may create a convexity condition where one does not 


exist without the terms. 
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b. If for a fixed Cy given by equation (3.9), the 
expansion (3.10) of a constrained state or control is inad- 
missible by a finite amount a. > 0) at telt T), its 
associated penalty functional (3.6) is unbounded as K > ©, 
This means that as E poo.tuEccontrrPburpton of a penalty 
term (3.6) to the augmented performance measure for an 
inadmissible state or control becomes very large. Therefore, 
AI Ja is being minimized under the condition of ever increas- 
ing Rn the constrained states or controls must at least 
approach admissibility. 

c. If for a fixed e given by equation (3.9), the 
expansion (3.10) of a constrained state or control lies 
completely within the admissible region, its associated 
penalty funetional (3.6) has limit zero as K +o, The 
significance of this result is that penalty terms of the 
form given by equation (3.6) will add less and less to the 
augmented performance measure as E, is increased for states 
and controls that are completely admissible. 

2. Convexity 

Property a discussed above is shown here. The 
theorem to be proved follows. 

Theorem 1. If a constrained state or control y(t) 


is bounded for t e Lt, »T] and is given by 


y (T)-y(t,) M mn(t-t_) 
y(t) en) sns cQ * 248 sn (3.10) 


where £s € Ro is defined by 


36 : 





c t [a a 
^y T Eus 


) 8» Y(t,), y(T)] £0 , 


(3.9) 
then the penalty functional 
Trey(t)-y,-y, ]2K 
J (y,K ) = m f M D (3.6) 
O 


n 
is convex on R. 


The constants Y y and Yr, define the 
admissible region for y(t) and r is a constant. 


Proof. Consider the case of Ky = 1. Equation (3.6) 
becomes 
A s S yc yes) dt. (3.11) 

Let 

A Y bay 

ai E (3.12) 
and 
A 
Hie (3.13) 


Substituting equations (3.12) and (3.13) into equation (3.11), 
we obtain 
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Cy 
| 


T 
2 
TEC J ly(t) - al? at (3.11) 
O 


ll 


T 
ro S [y (e) - ay(t)a + a] at . (3.15) 
O 


Substituting the expansion (3.10) into equation (3.15), we 
obtain 


y(T)-y(t,) mm(t—t,) )32 


"m 


(3.16) 


| y(T)-y(t_) M um(t-t_) : 
ay) t ELO (t-t ) t Ze am sin = | d + ‘| db, 


Equation (3.16) may be written as a quadratic functional 


of the form 


T 
doar, f [E4ej89) 008,5 * <ey.a(t)> + slat (3.27) 


ES 


where AER is 





T n(t-t,) sinen(t-t ) sinMn(t-t ) e tt 
of t ) = -2a|sin em es > — S Joe sm. 
S It, ? T DA T to T t. T t. 


(3.18) 
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Q (t) is the outer product given by 


T 
tem alt) (3.19) 
ha 
and 
A 
8- a^. (3.20) 


At an arbitrary fixed time telt oT], the integrand of 


equation (3.14) is 
Ly(ty)+al’= 3 CeyQ (tale >+<eyalty)> + 8+ (3.21) 


The first term on the right side of equation (3.21) is 


I y(T)-y(t,) mm (tt) 
2 eS 60€ = KS cot eee b DP | 


(3.22) 


Since the terms in the finite expansion of y(t,) given in 
equation (3.10) are linearly independent and analytic, y(t.) 


is different from zero almost everywhere and 


x : 
y (t) >0, ty € [t,,T], oc, # 0. (3.23) 
Thus, 
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and Q, (ty) is, therefore, positive semi-definite, at least, 
and positive definite almost everywhere for any E 7 0. 
Applying Theorem 4.5 of Reference 17 (p. AA the function 
given in equation (3.21) is convex on R" at t - Ec 

Next, consider the case where E is any positive 
integer. In Appendix E the following theorem is proved: 
if f(x) is convex on R” where X € R and f(x) 2 0, then 
f(x) is convex on R”? where K is any positive integer. 
Since 


2 


[y(t4) - a]^ » 0, (3:25) 


is convex, it follows immediately that 
| 2K 
[y(t,) - a] P (3.26) 


is convex on R” at a fixed time ty E [t ,T]. Since y(t) is 


bounded by assumption for t € [t >T], it follows that for 


llet f be a twice continuously differentiable real- 
valued function on an open convex set c in R . Then f is 
convex on c if and only if its Hessian matrix 


2 
_ EX un 
Q. E (24 400), Gs 3 6x) = 96,95, (Ei, eee 3 En? 


is positive semi-definite for every x € c. A quadratic function 


DO 3 Cx. Ox A ay to 


where Q is a symmetric n x n matrix, is convex on R” if and 
only if Q is positive semi-definite. 
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any finite positive integer Eos the expression (3.26) is 
bounded for t, € [t ,T]. By Theorem 4, (p. 536) of 


Re ference 18? 


T 2K 
a D T MEE (3.27) 
O 


is convex on R, Hence equation (3.6) is convex on RP for 


all finite positive integer values of Rn 


Bee Behavvor obetre-sNew-»Penalty Funetional*for an 
Inadmissible Constrained State or Control 


Property b is shown below. 


Theorem 2. Assume y(t) is bounded and given by 
equation (3.10) fort e [t T] where po #0, as defined by 


equation (3.9). If for a given y(t) 


y(t4) » yy * €, (3.28) 
or 
y(ty) & yy - €, (3.29) 
thet 
I p(x) = x) dt. 


Let T be of finite measure. Let f(t,x) be a finite convex 
function of x for each t and a bounded measurable function 
pae teren each x. Then T„ is a well-defined finite convex 
faction on L (T) which fs everywhere continuous with 
respect to thë uniform norm. 
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at some time t, e (tT), where ss: 0, then 


Ec = 0 B 
an p» Ky) (5 0) 
p 
where 
J E = — P at. .6 
NEN a) | E | (3.6) 


In order to prove this theorem the following lemma is 


required. 


Lemma 1. Assume y(t) is bounded by 
-> < M; < y(t) ul, = e c (3.31) 


and is given by equation (3.10) for t e [t >T] where 


ES X 0 as defined by equation (3.9). Then, if 
y (i) > Yu tT € (3.28) 


p 


at some time t, c (t T) for any En > 0, there exists a 


6 > 0 such that 
Fp 
y(t) >¥y +7 (3.32) 
for the finite time interval 


pU Era LEE E (3738) 


JB— < - o 


— M 








Similarly, if 
y(t) <y - € (3.29) 


at some time t, ce (t ,T) for any D 0, there exists a 


ô > 0 such that 
BD 
for the finite time interval 
Cy —- Š < t < ty t ó. (3.33) 


Proof of the lemma. First, it is necessary to show 
that the coefficients in equation (3.10) are bounded; that 


is 


| a, | m=] o SPM (3.35) 


= Mos 


where M3 > 0. To this end consider 


2 T y(T)-y(t ) mn(t-t. )32 
Jf on dt = y pA t-t ME a ¡Pin | dt 


O O 
(3599) 


( 6,5858,» (3.37) 
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where = is given by the definition (3.9). Q. is a 


symmetric matrix given by 























T- m = 
Co Ope. WIEN 2(T t0? 207 t.) 
2 T T 
0 Dt wm 2(T-6 ) _ a(T-t_) 
2 2T 2T 
Lo e e. e 4 e 
0 0 EE u) + 2 (Tae) 
2 MT MT 
s “Gr I 2(T-t,) DE DM 
T aT MT 3 3 
eU MO ETE E Dope JE 
T 2T Mn 3 33 
(3.30) 


where the terms with * are positive if M is odd and negative 
if M is even. Since y(t) is bounded by assumption and the 
expansion (3.10) is the sum of M * 2 linearly independent 


terms, we have 
5 2 
0 < [y(t)J dt < M SON) 
J 


where M, > 0. Using equation (3.37) and inequality (3.39), 


we obtain 


0 < (6,5858? < My. (3.40) 
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Q. is, therefore, positive definite. Using Theorem 2.5 


of Reference 15 (p. beo. we have 
(e,,Q,o0 9 » X. |le E (3.41) 
AY” Y 


where A > 0 is the smallest eigenvalue of Q, and 


Ile i] = V< Ey ey . Therefore, using the inequality 
(3.40), we obtain 


M 
2 h 
Hell s. (3.42) 
Since an? m= 1,2,...,M is a subset of gs it follows that 


la. | «M m=1,2,...,M (3.35) 


“al 3 


where Ma 2 0. 


Now consider the derivative of equation (3.10) 


which is 
; y(T)-y(t_) M a mm(t-t ) 
y(t) = ——— E Ar CS e (3.43) 
T-t_ da m I-t, T t. 


liet Q = RITE be a symmetric n x n matrix. Then Q is 


positive definite if and only if there is a k > 0 such that 
2 
<v,Qv> > k]|v]] 
for all v in R^, where ||v]] = ,[(V,V> is the Euclidean 
Horm of v. 
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Taking the absolute value of both sides of equation (3.43) 
and applying inequality laws, we obtain 
y(T)-y(t_) M nn(t-t_) 


ly(t)| < | ae ] + | Da, 9 cos = S| (3.44) 
=o NICE D 


which further simplifies to 


: y(T)-y(t_) M 
KEDI < rr > lanl r= - (3.45) 
O m=1 O 





Applying the inequality (3.35), we obtain 


Ber. y(T)-y(t_) MoT M 
nce) |< — «s T 2,m (3.46) 





M (3.47) 


|^ 


for allt e (to, T) where Mo > 0. The first part of the 
lemma as expressed by the inequality (3.32) will now be 


shown. Let 


s P ` (3.48) 


as shown in Figure 4. Consider 


t. 
y(t) = y(t,) + f Code (3.49) 
£ 


* 
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Figure lj 


Constrained Variable vs. time 
Applying inequality (3.47), we have 
y(t) 2 y(t,) = Me [ts x t | . £515 519) 


Consider t in the interval t, - ô <t < tą t ô. In this 


interval ó satisfies 
6 > |, — t|» (3.51) 
Applying inequality (3.50), we have 





Using the definition (3.48), we obtain 


£ 
y(t€)»y(t)--. (3.53) 


Applying inequality (3.28), we obtain 


- > 
AER e (3.54) 
or 
€ 
y) > yu t — ° CE 


thus proving the first portion of Lemma 1. The second portion 
of the lemma as given by inequality (3.34) can be proved in 
a similar manner. It is possible at this point to return 


to the proof of Theorem 2. 


Proof of Theorem 2. If inequality (3.28) applies, 
then by Lemma 1 inequality (3.32) is true. Considering the 


integrand of equation (3.6), since Yy ? Yi? it follows that 


, € 
ie Pj y - 
= Y y E : EZ Y y > ae 


YM YI YM YI 


Irt + 6 = t < t, t ô. Purther simplification yields 
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zer: 


€ PY MI, ae 
YMYL 


YMIL 


|v 


€ 32K 
«gi 
YML 


(3 


(3. 


Since the integrand of equation (3.6) is nonnegative for 


tt (t T) and r > 0, it follows that 


rf E dt 
6. YM YL 


K 
J (y, p! 


+6 

tx 2y (t)-y y, 2K 
> e f | ee Cum P at 
"d rs In Yi 


Applying inequality (3.57), we have 


t,tó 
INY K) >r J E on. dt 
tg-6 YM IL 


The integration of inequality (3.59) yields 


i € 
(y K > r È +e] P 28 
p p A Yu L 


Since yu > Yk and en 220x115 Lollowswenat 


II 
8 


Limit 1 4+ Ep eK, 
E YMYL 
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58) 
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and in view of inequality (3.60) 


Limit zuge 
i Jy fy sk.) a ° (3.62) 


If inequality (3.29) applies, then by Lemma 1, 


inequality (3.34) is true. Since Ym > Yi» 1t follows that 


€ 
> < E e 


= = (3.63) 
mn cm 
Tor Ty - 6 < t < ty t ó. Simplifying, we obtain 
2y(t)—-yu-y € 
(A i | | (3.64) 
NEU TEA 
Squaring both sides òf inequality (3.64), we have 
2y(t)-yy-y, 77 e. 72 
| J Es > | . (3.65) 
IM ’L YM YL 


Raising inequality (3.65) to the K. th power, we obtain 


2y(t)-yy,-y. 12K € 2K 
— Pr [1 2] P (3.57) 
M "L ML 


which is identical to inequality (3.57). The remainder of 
the proof follows inequality (3.57) to equation (3.62) 

exactly. The theorem is proved for the open interval te(to T). 
The extension of the theorem to cover the closed interval 


L£ BET] is not difficult: 





H. Behavior of the New Penalty Functional for an 
Admissible Constrained State or Control 


Property c above is shown below. 
Theorem 3. Assume y(t) is given by equation (3.10) 


forte Lt TJ. If for a given y(t) 


yy te, < y(t) « yy - Eq (3.66) 


for all t e [toT] where a 0, then 


Limit 


By > 


SO .6 
JS D. (335670) 
Proof. From inequality (3.66) it can be seen that 


yr + ek (3.606) 


Since yy » yj, inequality (3.68) may be rarranged to the form 


2€ 
je E O) 


YMYL = 


3 (3.69) 
The inequality (3.69) will become useful shortly. 
Starting with inequality (3.66), multiplying through 


by 2, and subtracting Ym and yp» We obtain 


2(yp te )-Yy7 YT, < 2y(t)-yv Yr $ 2 (Yy Eg) IMI ° (3.70) 


al 








Dividing through by the positive quantity Yu T Yg» We obtain 


2(yr*e )-Yy7 Jr, 
Eu ov « 


YMYL 


Iu YL 


2y (t)-yu- Yr, 2(yu-e )-yy YL 


< —— L a, 
_ Iu L 


This inequality can be reduced to 


2€ 
ae oa 
Iu YL 


In view of inequality (3.69), inequality (3.72) may be 


rewritten as 


2y(t)-yy-y DE 
|: i 
= JM YL 


y A 
YMYL 


2y(t)-y,-y 2€ 
I eim F 
YM Jr, e JM Jr, 





Raising inequality (3.73) to the 2K power, we have 


|^ 


Frl 
MIL 


Observing that 


A 
E 


we have from inequality 


Eo y 
YM YL E 


for all te [toT]. 


easily seen that 


2€ 2K 
2 - | p 
YM YL 


= ° ° 
= 3 
YMYL 


(3.78) 


Ë E 2€ p 


YM YL 


In view of inequality (3.76), it 


Dc 


(3.71) 


(Berne) 


(3.73) 


(3.73) 


Nc 


SE 








T T 

2y(t)-y.,-Y- 22K Qe 2K 

f |——— " P at < J 2 ee) 
c. YMYL t L 


Performing the integration, we have 


T 
2ay(t)-y.,-y. ]2K 2€ 2K 
f E P ac < Ë e | P (TRER aoa OD 
Ur YMYL YMYL 


In view of inequality (3.69), and the fact that Yy ? Yp 


it follows that 


E 2K 
uut. eu (3.79 
Rp ^^ ea | 


Observing inequality (3.78), we obtain 


s | 

2y (t)-y.-y. ]2K 

Limit a p -- (3.80) 
| e Jw TL 


The theorem is proved. 
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IV. THE ALGORITHM 


This section describes the algorithm used for minimizing 


the augmented performance measure. 


A. GENERAL MINIMIZATION STRATEGY 
1. Sequence of Unconstrained Sub-problems 

A sequence of unconstrained sub-problems is solved 
by the algorithm. In each sub-problem the algorithm mini- 
mizes the augmented performance measure for given values 
of the weighting factors (e and r) and the inequality con- 
straint penalty term power (K.). After an appropriate 
stopping criterion is satisfied, e is reduced, Rn is 
increased, and a new sub-problem is commenced using the 
optimal solution to the last sub-problem as a first guess. 
This procedure is repeated until enough sub-problems are 
completed to meet a second stopping criterion. 

The algorithm is programmed to do one sub-problem 
on each computer run. The results are stored on an external 
storage ev between computer runs and are retrieved at 
the commencement of the next run (new sub-problem). 

2. Minimization Strategy 

The algorithm minimizes by either the FNR or MNR 
method. The user must decide which method to use on each 
iteration. This is a matter of experimentation, especially 
for the first two or three sub-problems. An effective pro- 


cedure is to run the sub-problem once using the MNR equation 
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throughout and once using the FNR equation throughout. From 
these results an effective minimization strategy can generally 
be deduced for the sub-problem. Occasionally further exper- 
imentation is required. This experimentation points out 

the advantages of using separate computer runs for each 
sub-problem. Once a sub-problem is completed and the results 
stored, the computation does not have to be redone each 


time an experimental run is made in the next sub-problem. 


B. COMMENCING THE PROBLEM 
l. Initial Decisions 

Three interrelated decisions must be made to begin 
a problem. First, the number of time points K must be chosen. 
Second, the number of coefficients M for each state and 
control expansion must be chosen. The same number of coeffi- 
cients is used for all expansions in a given problem in 
this dissertation, but this is not a requirement. From a 
theoretical standpoint it is desirable to use a large number 
of coefficients and time points to insure that an adequate 
approximation of the optimal control and state trajectory 
iS Obtained, but practically, computer time and storage 
requirements Limit the number of each. The computational 
penalty for using a large number of coefficients is the 
more severe of the two as the number of equations in (2. 36) 
and (2.37) which must be solved is equal to the total number 
of coefficients plus the number of free end conditions. The 
solution of equation (2.36) or (2.37) represents a considerable 


portion of the overall computer time. 
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The third decision involves the initial values of 
€e; r, and oo The weighting factor r for the inequality 
constraint penalty terms is held constant throughout the 
entire problem. A satisfactory value used in all problems 
in this dissertation for all inequality constraint penalty 
terms is r = 100. An initial value of Ls which has worked 
well in all problems is K. = 4, Larger values of K. gener- 
ally cause computer overflow in the first sub-problem. 
With these values chosen there exists a region of e's for 
which the first sub-problem will respond to an appropriate 
minimization strategy. This acceptable range of starting 


e's is different for each problem but is in the range 
10, Da Io 


for all problems solved herein. Numerical experimentation 
is the only method available to determine an acceptable 
starting e. There is no theoretical requirement to use 
the same value of e for each state equation equality con- 
straint term in the augmented performance measure or the 
same K_ in each inequality constraint term, but the use of 
different e's and K's has never been required. 
2. Initial Guess for the Unknowns 

Once the above three initial decisions are made, an 

initial guess for the vector of unknowns c is required. The 


vector c includes all coefficients and free end conditions. 








All coefficients are set equal to zero initially unless 


there is good reason to make a different choice. 


C. ITERATION 
l. Required Vectors and Matrices 

The states and controls are calculated at each time 
increment by evaluating the functional expansions, The W 
vector defined in (2,15) is calculated using these states 
and controls, Next, the gradient matrix (2,34), the aug- 
mented performance measure (2,16), the symmetric matrix 
(Vow) A a aná the vector (VW). m, are calculated. 
E At this point the algorithm begins the iteration 
process with either the MNR or the FNR method depending on 
the value of a flag set by the user (the method selected 
is based on the iteration number being performed). If the 
MNR method is to be used, equation (2.37) is formed. If 
the FNR method is called for, the three-dimensional array 


(2.35) is calculated and the symmetric matrix (V. “We M 


^ 


is formed. It is prohibitive to store the Eu ee secl 
dimensional array, but a feasible alternative is to multiply 
each matrix in this array by was the matrix is calculated 
and store the Eee column vector, Once a matrix in 

the three-dimensional array is multiplied by w, it is no 
longer required by the algorithm. The next matrix in the 
array is calculated and stored in the same storage locations 
used E the first matrix. Only the symmetric matrix 


(V. ^w), Tw. need be stored. The total increase in storage 


^ — ^ 





requirements of the FNR method over the MNR method using 
this computation technique is less than 10 percent in the 
problems solved herein, It is also imperative in terms of 
computation time to take full advantage of the symmetry of 
the matrix (V Aw) Tw. Due to this symmetry it is necessary 
to calculate BEA ome column of the first matrix in the 
array, two columns of the second matrix, and n columns of 
the neh matrix. By taking advantage of the symmetry the 
average time for each FNR iteration is approximately twice 
the time for each MNR iteration. 
2. Solving the Linear System 
At this point equation (2.36) or (2.37) is formed 


and must be solved for Ac. This is a linear system of the 


form 


t >> 


= (115 


and is solved in the algorithm by one of three methods 
available to the user. They are: 

a. Gauss elimination with improvement by residuals 
using total pivoting, | 

b. Gauss elimination with improvment by residuals 
using main diagonal pivoting and a computation technique 
which capitalizes on the symmetry of A, anda 

c. Gauss-Seidel iteration. 

In the problems solved herein the number of unknowns 


varied from 37 to 74. In spite of the large number of 
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unknowns involved, the elimination methods required less 
computation time to solve the linear system than the Gauss- 
Seidel iteration method. It was observed for the problems 
solved that total pivoting was not required in the elimina- 
tion method. Method b, therefore, was the most economical 
and effective method for solving the linear system and was 
used in all problems, Method c is retained in the event 
that the algorithm is used to solve problems with a larger 
number of unknowns. 

In each solution of equation (4.1) one improvement 
is made using residuals. That is, after equation (4.1) is 


solved, 


¿> 


x-b=r (4.2) 
is formed. The system 
Ay =m (4.3) 


is solved and the resulting y is subtracted from x to form 


the final solution to equation (4.1). 


B Interpolation 


Tabular functions of two independent variables are 
used extensively in the problems to represent aircraft and 
missile parameters accurately. Parabolic interpolation is 
used to obtain the functional values in these tables and 


the required first and second partial derivatives. The 
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derivation of the necessary difference equations for para- 
bolic interpolation in two independent variables is presented 
in Appendix C. 

Excerpts from the tabular data used in the problems 
is presented in Appendix B along with graphical representa- 
tions of the data. The data represents typical supersonic 
aircraft and missile performance parameters and has been 
obtained from several sources. Considerable effort was 
expended to smooth the data before the tables were constructed 
since finite difference methods were used not only for 
funetional values but also for first and second partial 
derivatives. 

H, Stopping Criteria 
Once the linear system is solved, a new c vector is 


calculated from 


Il i i (4,4) 


20 


At this point a stopping criterion is tested. If 


i 


IT, 


ER, (4,5) 

the sub-problem is finished. Otherwise the iteration process 

is continued. At the completion of the sub-problem the 

results are stored off line. The computer run is complete. 
To begin a new sub-problem a new computer run is 


initiated, recalling the results stored from the last 
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sub-problem, Epsilon is decreased, Eo is increased, and 
the minimization strategy is altered by the user as required. 
Typically, g is divided by a factor of between two and ten 


and 2 NosNUunePCased" by "Cwo"OP"fouPT" "T58tC. is, 
K = K ` + orj. (4,6) 


More ambitious policies usually result in failure of the 
algorithm. 

Sub-problems are solved until a second stopping 
criterion is satisfied. Several criteria are possible to 
end the problem. A method used successfully involves 


observing 
JI. +d (4.7) 


and stopping when this sum, which represents the penalty 
terms due to the equality and inequality constraints without 
weighting factors, ceases to decrease significantly between 
sub-problems. 
5. Flow Chart 
A flow chart of the algorithm is given in Figure 5. 
6. Integration 
At the completion of the last sub-problem a check 
on the degree.of satisfaction of the state equations is 
obtained by comparing the state expansions with the state 
trajectory obtained by integrating the state equations with 


the control expansions. 


61 


MI ¿A o 
> © q o 
— - > 


d D uud 2— S 





Read constants, flags, initial guess for c, initial conditions, 
weighting factors, inequality constraint power, 
fixed final conditions, stopping criteria. 


Evaluate functional expansions at each time point 


Calculate w, Vow 











Form Yan (vw) w, and A = (Vw) 


FNR Calculate N W 













Solve for Ac = - A (Vw) W 


i+] i 
Calculate new c vector, c = G + Ac 


Yes 


Sub-problem done. Reduce e's and increase Kb S 


o E E (9 0%] < 52? = 





Figure 5 
Algorithm Flow Chart 
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V. A MISSILE INTERCEPT PROBLEM 


In this section a short-range missile intercept problem 
is solved. An air-to-air missile launched from an attacking 
airplane is to intercept a constant-velocity target in 
minimum time. The missile is restricted to move in a plane. 
The orientation of this plane is defined in three- 
dimensional space as the plane containing the position of 
the missile at launch, the position of the target at launch, 
and the velocity vector of the target. The assumptions 
applied to the problem, the coorinate systems used, the 
nomenclature, and the derivation of the equations of motion 


are presented in Appendix A. 


A. PROBLEM FORMULATION 


1. State Equations 


The state equations derived in Appendix A are 











° gTh gpsa 5 gosa 5 I gW. l 
M = aW Th cosa - ^W M C cosa — 5; M Cysina S sind 
e e e e 
(5.1) 
° eThy Th sina gosa gpSa gW. cos9 
Lu m 2 zw. Mc,sina * sy N o> av, a 
e .e e 
(522) 
X = aMcos9 - aM ncosYy (5.3) 
Y= aMsind - aM, siny ; (5.4) 
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The states are Mach number M, flight path angle 0, relative 
range X, and relative cross-range Y as defined in Figure 40. 
The control is angle of attack q. 

The normalized missile thrust Th is considered 
constant until missile engine burnout t 


B and zero thereafter. 


That is, 


Th 


l 
= 


5 te [0,tg) (5.5) 


Th = O st E (t5, T] . (546) 


Other parameters considered constant are 


32.1725  ft./sec.^ 


g = 
Thy = 3500 TDS; 
TN 
W = 200 lbs. 
e 

= 0.35 ft.* 
t5 = 7 sec. 
a = 1077.0 ft sec: 

0.001756 slugs/ft.? 


D 
tl 


The values of the constants in equations (5.7) are based on 
an engagement at 10,000 feet altitude. It is convenient to 


group these constants as 


| etn | 
M 
a (5.8) 
d aW, 
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In addition, 





- (5.9) 
Ma 

— B. 
aw (5.10) 


the computer solution is aided by normalizing 


the states X and Y by defining 


II > 


I> 


= (5.11) 


Y 
m (5.12) 


where X and Y are assigned nominal values of 10,000 feet. 


pee 
II 


are 
and 


air 


E 


= 


><] 


<| 


2; 


With these adjustments the state equations ane 


Thcoso - 


Th sino 


M 


M cos - 


M sino - 


2 2 : 
CM C Cosa - CM Cu Sina - C.W sine (5.13) 
W,cos9 
- CMC sina + C¿MC, cosa - C. — (5.14) 
Š M_cos (5.15) 
x Ten 
2 Masin (5.16) 
y Mm Y : s 


Tabular Functions 


The axial and normal force coefficients Ca and Cn 


given in Appendix B as tabular functions of Mach number 


angle of attack. This data is based on a typical air-to- 


missile. 
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3. Performance Measure 


The performance measure for this problem is 


T 
J = T dt 
0 


H. Inequality Constraints 


Four inequality constraints are required. 


of attack must satisfy 


I 
SE 


< Q < 


ja 


(5.17) 


The angle 


(5.18) 


From structural considerations the load factor must satisfy 


ET 


5S: End Conditions 


(5.19) 


In order to describe the initial conditions for the 


problem it is necessary first to pose the problem in the 


(X,Y,2) coordinate system shown in Figures 41 and 42 of 


Appendix A. The problem chosen for presentation in this 


section involves a target in a shallow climb at short range 


with a slight altitude disadvantage crossing the attacker's 


flight path extension at 90°; i.e. 


Rn = 15000 ft. 


Am = -3000 ft. 
ad ° 

Bm = 90 
= ° 

Om 10 
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(5.20) 
(52) 
(T22) 
(5123) 








Following the procedure outlined in Appendix A, the 


remaining unknown parameters are 


Ma 7 0.8 (5.24) 
= 42,35° (5.25) 
a 51.526 lbs. (5.26) 


The initial conditions as computed by this procedure are 


M(0) = 0.8 (5.27) 
8(0) = -49.573° (5.28) 
x(0) = O (5.29) 
peo) = 0 a (5.30) 


The final conditions as computed by this procedure are 


0.9920 (ora 315) 
EJUS!" . (5292) 


x(T) 


y(T) 


Note that x and y are the normalized relative range and 
relative cross-range, respectively, as defined by equations 
(5.11) and (5.12). The states X and Y in equations (5.11) 
and (5.12) are defined as the position of the missile in 
the (X,Y) coordinate system shown in Figure 40. At t = 0 
the missile is at the origin of the (X,Y) system, hence 
equations (5.29) and (5.30) apply. At t = T the missile 
must intercept the target. Since the (X,Y) system becomes 


fixed with respect to the target at launch, x(T) and y(T) 








are equal to the normalized coordinates of the target in the 
(X,Y) system at launch and are given by equations (A.93) and 
(A.95). 


B. THE EPSILON METHOD FORMULATION 
1. The Augmented Performance Measure 
Using the penalty functionals described in Section 
III for inequality constraints, the augmented performance 


mr 


measure is 





T 
+ 2 f [à - o mos  C.M*C. cos + CM sin + oM sino] at 
e+ 1 A a E OS 
| 


š W cos0 42 
1 f Thsina we c ] 
+ 2 d ð -C) M + CMC psina CMC cosa t C3 M dt 


T 2 
° a a ' 
E - — Mcos8 + — cosy | dt 
Xm ES) 
E 2 
È A y ine + 9 Mysiny] dt 
0 





l 2K Tr 2, 42K 
ee | p f E | p 
+ rf |a SUE ie i 5088 dt 


where € and r are weighting factors and ô is a constant used 
to make minor adjustments in the admissible regions. In all 
problems 6 is given a value of 1.03. This adjustment is 


applied to all admissible regions of constrained states and 
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s 





controls. The power ds is limited in computation to approxi- 
mately 30 before computer exponent overflow problems develop. 
With this value of Ko it is desirable to adjust all admis- 


sible regions slightly as can be seen from Figure 6. 


E ius c. 60 
d Yu YL 
| 6 YM YL | 
| | 
i 2 
| | 
A| | 
K_=30 | K =30 | K =30 
621.03 jl ô=] | 6721.03 
| y 
| 
| 
Y YM 
Figure 6 


Adjusted admissible regions 
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The required elements of w are 


M 2 2 
w. = IM, C,Theosa, * CM. C, cosa, us Cy Sina, 
(5.34) 
11 
+ CW sine, | |£] a Pisa, K 
" Thsino, 
Week = lo, = Ci ME + d eU: - BE ili 
(5.35) 


W coso 1 
y + Y + | | A= 9 k = Dee ies. aig K 


1 
a “4 


E- Ws a At = 
Wont = | - 7 M,cosé, +7 M„eosy| [| A. E COSO) 


` i . 1 ; 
Wokx+x = Ë - 5 M, sin0, +F Mpsiny] [22] k2, K CONES 


20, VK š | 
Wuk+k = = P (rat)  , ku D. ck (5.38) 
aĝ, M, JK 
= kk |p % 2 
WoK+k = ES (rAt) , k = A S (5.39) 
% | h 
Worst = (Ck dA (5.40) 
where 
M. A M(t) - eto, 
t=(k-1)At 


TO 








2, Functional Expansions 


The state and control expansions written in discrete 

















form are 

M. = M; + Ka (k-1) + ^ a sin r y ce ES PI 
6, = 0, + A (k-1) * T b, sin nik) ee oe 
Een on Een 
Y, = y, + x (k-1) + E a sin cm M n ne 
a. = a, + Ey (k-1) + È e sin ss le a 


3. Vector of Unknowns 


The elements of the vector c are 


T 
ce = (84,85, «58g, b4 b5, «by 045055 «9C 


^ 


d415055* «50g 5915955 «ey Mgs eg tA At) , 


H. ‘Partial Derivatives of the Tabular Function 


The values of C, and C 
A N 5C 


by parabolic interpolation along with the values of 5 A 


< 
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en 
> Ya 


(5.41) 


(Se) 


(AS 


Car 


(5.45) 


(5.46) 


are obtained from the tables 


dC 


3 














C aC, EU O or ace 
— , — and 
3M^ jae” Ma 3M XE tete ° 
ge 
SMSa ` The procedure is outlined in Appendix C. 


5. First Partial Derivatives 
The first partial derivatives indicated in equation 
(2,34) are required. These partial derivatives are easily 
obtained from equations (5.34) thru (5.40) by taking the 
partials of these expressions with respect to the vector c. 


The expressions are too numerous to include. A typical term 


is 
2W, aw, aN Jw, 9M 
et (5.47) 
m oMx m k m 


For 1 < k < K, w, is given by equation (5.34) and the partial 


k 
derivatives indicated above are 


JW | 
E- [=] (5.18) 
oM,- € 
OW), > me 
3M. = EACUS + C.M, EN cosa, + nla Ne ale 
d Ch | (5.49) 
2 k At 1% 
+ C M, s= sina, | | E | 
3M 
k _ mm mm (k-1) 
ga, — (K-1)At KET (5.50) - 


T2 








— = sin MI . (5.51) 


Notice that C and C are functions of M, as well as a. 
A Ny K K 


6. Second Partial Derivatives 
The second partial derivatives indicated in the three- 
dimensional array (2.35) are also required, Again the expres- 


sions are too numerous to list. A typical term for 1 < k < K 


is 

3^w ow 

Ja k - A EJ NODE 

gm x m 
ow, A 
Letting coc R „were Mave 
m 

3^w 3M oM 
ek OR -3R K BRO CL (5.53) 
92,08, ga, OM, da, 9M, 98, 





oR = 0 (5.54) 
ð Mk 
3C ale 
pw o.c. co + 4CM iio o + O.M ° “k cos 
aM 7 ae Sr MOMO Ser 2" k 2 ak 
k k k aM, 
aC 950 
Ny > N, 
+ ¿Con sina, + a CAM, aM sina, t C M, Fr sing, | 
k k 3M 
K 
[At] tsin zk (5.55) 
€ K-1 


=> \ 











@ 
N 


_ £T £T(k-1) 
da, ^ (XC-1)8t 799 "X-1 (5.56) 
oM 
i sin Enel : (Sa 
C. RESULTS 


The problem was solved twice: once using 8 coefficients 
for each expansion (problem A) and once using 12 coefficients 
for each expansion (problem B). In both cases K = 21 time 
points (20 time intervals) were used. The initial guess 


for the ec vector was: 1 


all expansion coefficients = O 
HE = 1,4 
6(T) = 0° 
a(0) - 209 (5.58) 
a(T) = 0° 
At = 7/20 sec. (T= 7 sec.) 


l. Problem A - 8 Coefficients for each Expansion 


Four sub-problems were required to solve the problem. 
Table 3 gives the weighting factor values, optimization 
strategy, and computer time for each sub-problem. Table 4 
gives the components of the minimum augmented performance 


measure for each sub-problem where 


aa =J + rJ (5.59) 


fe 


sub- € K I* | optimization Copal 
problem P strategy** time 


io on MMMMMMMM 2128" 


0.67 x 1072 FFFFFF 3! 38" 
DS x JG ° FF 154" 
0.5 x 107? F 1'32" 





* number of iterations required 
** M - MNR method 
F - FNR method 


Table 3 


Weighting factors, optimization strategy, and C.P.U. time for 
missile intercept problem A. 


Sub- 
problem 


4 
4 
4 
4 


0.6853 x 10. 0.2551 


0.4626 x 10° 0.1287 
0.4624 x 10° 0.4619 


0.4624 x 10° 


0.3691 





Table 4 


Components of the minimum augmented performance measure for 
missile intercept problem A. 
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Figure 7 is a plot of the augmented performance 


measure vs, iteration number, 


Augmented performance measure 


A o 


SUB-PROBLEM 


0-0 Og 66 B-0-0-0 
2 1—4 2 —k3 44| 
SUB-PROBLEM 





0 2 4 6 8 10 12 14 16 18 20 22 24 
Iteration Number 


Figure 7 


Augmented performance measure vs. iteration 
number for missile intercept problem A 
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Iterations performed by the FNR method are indicated by a 
solid line. Iterations performed by the MNR equations are 
indicated by a broken line. The figure shows several signi- 
ficant characteristics: 

a. the failure of the FNR method on the first itera- 
tion of the problem (the initial guess is too far from 
optimum to allow the FNR method to converge); 

b. the superior terminal convergence produced by 
the -FNR method close to the minimum; 

c, the oscillatory results produced by the MNR 
method as the minimum is approached. After the commencement 
of sub-problem 2 as shown in Figure 7, the FNR method is 
used exclusively to the end of the problem. For comparison, 
sub-problem 2 is commenced with the MNR method and allowed 
to run to 23 iterations. At the beginning of sub-problem 3, 
Ja increases slightly and never returns to the minimum value 
obtained in sub-problem 2. This is an indication that e has 
reached a point where smaller values have little influence 
en the reduction of Joe 

Figure 8 is a plot of the performance measure vs. 
iteration number. The performance measure (final time) 
increases with iteration number. As the algorithm iterates, 
At is being increased to reduce J. (the term in the augmented 
performance measure reflecting the degree of non-satisfaction 
of the state equations) while insuring that constrained 
states and controls are kept within admissible bounds by 


Beducing or holding down In‘ With small values of e and 


Tr i 








Sub-problem number 
7.8 Í 2 34 





A 
o. 


` 
+ 


Performance measure - T(sec.) 


0 2 4 6 8 10 12 1 16 18 20 
Iteration number 


Figure 8 


Performance measure vs. iteration number 
for missile intercept problem A. 


large values of Ko the end result of minimizing the augmented 
performance measure is a control history and state trajectory 
which minimizes the time to intercept while satisfying the 
state equations and inequality constraints; all to a 


reasonable degree of accuracy. 
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Figure 9 is a plot of Jo and J. vs. iteration 


number, 


Sub-problem number 
1 2 3 4 





0 2 4 6 8 10 12 14 16 18 
Iteration number 


Figure 9 


J. and J vs. iteration number 
for missile intercept problem A. 


After K is increased at the beginning of sub-problem 2, the 
inequality constraint penalty terms J become very small. 
This is because the angle of attack is well within its 


admissible region. 





It is evident from Figures 7, 8, and 9 that two 
sub-problems are sufficient to obtain a reasonable solution. 


The overal stopping criterion 


“y 1 *,1+1 6 (5.60) 


| w: Be Ios 
S p S p 
where i is the sub-problem number is satisfied after the 
third sub-problem. However, at this point the value of K. 
is 8 which is not large enough to provide desirable 
penalty functionals for the inequality constraints. It is 
necessary to provide as large a value of K. as the computer 
will allow for the final sub-problem to insure that the 
minimization is not influenced by the inequality constraint 
penalty terms when the constrained states and controls are 
within their admissible regions. Accordingly, a final sub- 
problem 15 performed with K. = 32, The algorithm is able 
to handle the increase in K from 8 to 32 in one step only 
because no constraint boundaries are active in the solution. 

Figure 10 is a plot of the angle of attack expansion 
computed at the end of the last sub-problem. The region of 
admissible angles of attack is shown. 

Figures 11 and 12 are plots of Mach number and 
flight path ance. time. In each plot two curves are 
shown; one is the expansion for the state as computed at the 
end of the last sub-problem; the other is the state trajec- 
tory obtained by numerically integrating the state equations 


with the optimal control expansion. 
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Figure 10 


Angle of attack vs. time for 
missile intercept problem A. 
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Mach number 


integration 


expansion 





i 2 3 4 5 6 7 
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Figure 11 


Mach number vs. time for 
missile intercept problem A, 
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Flight path angle (degrees) 


Time (sec.) 


integration 


expansion 





Figure 12 


Flight path angle vs. time for 
missile intercept problem A. 
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Figure 13 is a plot of relative range X vs. relative 
eross-range Y. As before two curves are presented: one 
represents the expansions of the states as computed at the 
end of the last sub-problem; the other is the state trajec- 
tory obtained by numerically integrating the state equations 


with the optimal control expanslon, 


Relative range (ft. x 10°) 


0 2 á 6 8 10 12 


o 


4 
C 


integration 


Relative cross-range (ft. x 10%) 
3 


expansion 





-12 


Figure 13 


Relative range vs. relative cross-range 
for missile intercept problem A. 





Figure 14 is a plot of range X' vs. cross-range Y' 
where both quantities are obtained by transforming the 
expansions of the states X and Y obtained at the end of the 
last sub-problem from the (X,Y) coordinate system to the 
inertial coordinate system (X',Y') fixed at the missile 


launch point (Figure 40). 


missile 


Range (ft. x 10°) 





4 ó B 3 
Cross-range (ft. x 10") 


Figure 1! 


Range vs. cross-range for 
missile intercept problem A. 
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Figure 15 is a plot of load factor vs, time where 


the load factor is given by 


= 
enge (0 M) (5.61) 


and the states used in equation (5.61) are the state expan- 


sions obtained at the end of the last sub-problem. 


Time (sec.) 


Load factor 





Figure 15 


Load factor vs. time for 
missile intercept problem A. 
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Although load factor is not a state but a function of states, 
the plot is important because it shows that the load factor 
constraints as given by (5.19) are not active. 

It might be suspected that the optimal trajectory 
is a maximum performance turn limited by a constraint 
boundary (angle of attack or load factor) followed by a 
straight line path to intercept. This is not the case. The 
initial turn rate of the missile is small compared to its 
maximum turn rate capability. This is due to the high 
induced drag associated with high angle of attack turns 
which would reduce the missile's longitudinal acceleration 
capability. Also there is no straight line segment to the 
trajectory although the turn rate of the missile is very 
small at intercept. 

The first sub-problem was also solved with an 


initial guess for the c vector of: 


all expansion coefficients = 0 


M(T) = 2.5 l 
g(T) = 09 (2,825) 
al O0 = 108 
cU = OD 

At = 8/20 sec. (T = 8 sec.) 


The sub-problem reached a minimum of 7.544 which compared 
favorably with the minimum reached by the first initial 
guess given by equations (5.58). This gives an indication 


that the minimum attained is the global minimum. 
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2. Problem B - 12 Coefficients for each Expansion 





Four sub-problems were required to solve this 


problem also. Tables 5 and 6 summarize the performance of 


sub- r K I* opeimszstionlwC.,P.U, 
problem strategy** time 


1 1.0 x 107? [100 | 4 | 8 | MMMMMMMM 3' 31" 
2 0.67 x 102|100 | 6 | 2 | FF 
3 0.5 x 10? 1100 | 8 | 2 | FF 
h 0.5 x 107? 32 | 2 | FF 


2'20" 
19" 
2120" 
* number of iterations required 
** M - MNR method 
F - FNR method 


the algorithm. 






















Table 5 


Weighting factors, optimization strategy, and 
C.P.U. time for missile intercept problem B. 


022125 0.1849 x 10^ 


013565 u 0.3412 x 2 


0.1513 0.3474 x 1307 


0.1339 0.6139 x 1073 





Table 6 


Components of the minimum augmented measure for 
missile intercept problem B. 
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Figure 16 is a plot of the augmented performance 
measure vs, iteration number and corresponds to Figure 7 


for problem A. 
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Figure 16 


Augmented performance measure vs. iteration number 
for missile intercept problem B. 
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* 


A second failure mode of the MNR method close to the minimum 


is shown. The MNR method produces a divergent J. in the 


Second sub-problem. 


Figure 17 is a plot of range X' vs. cross-range Y' 


obtained in the same manner as in problem A (Figure 1H). 
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Figure 17 


Range vs. cross-range for 
missile intercept problem B. 
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A comparison of the tables and figures for problems 
A and B show that the optimal control and trajectory have 
not been markedly affected by the increase in the number of 
coefficients from 8 to 12 for each expansion. It is prohib- 
itive in terms of computation time and storage requirements 


to increase the number of coefficients further. 
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VI.. A CLIMB PERFORMANCE PROBLEM 


In this section a climb performance problem is solved. 
A supersonic fighter aircraft is to climb from sea level to 
high altitude in minimum time. 

Flight test experience has shown that to climb to 
altitudes on. the tropopause in minimum time a supersonic 
fighter must execute a maneuver which typically includes: 

a. a subsonic climb to an altitude near the 
tropopause; 

b. a level or near level acceleration to some 
supersonic Mach number; 

c. a "zoom" climb to the desired altitude trading 
kinetic energy for potential energy. This technique has 
been used extensively in the past decade for establishing 
climb records and in fighter-interceptor tactics to attain 
altitudes higher than the aircraft's service ceiling for 
short periods of time. 

Several factors contribute to the optimality of this 
type of maneuver. They are: | 

a. a fighter's maximum Mach number or "placard" 
limit which arises from dynamic pressure and/or thermal 
limitations and is a function of altitude; 

b. a fighter's transonie drag charaeteristios; 

c. air temperature variation with altitude; 

d. air density variation with altitude; 

e. turbojet engine maximum thrust variation with 
altitude, 

9e ' 








Optimization techniques were first applied successfully 
to this problem by Bryson [Refs. 1 and 2]. The method of 
Steepest descent was used successfully to predict the type 
maneuver described above for a typical supersonic aircraft. 

The epsilon method is applied to the problem herein to 
demonstrate the method's power, A direct comparison of 
methods is not made as the mathematical model used here has 


been.improved considerably over that used in Reference 2. 


A. PROBLEM FORMULATION 
l. State Equations 
The state equations for this problem derived in 


Appendix A are 








M = E cosg -= £ siny = 3j cM C, Cony) 
e 
Sa 
, . EIh sina SP, _ gceosy 
et ow ONL aM (6.2) 
h = 2 sin Y (6.3) 
Hi 


The states are Mach number M, vertical flight path angle y, 
and normalized altitude h, The control is angle of attack a. 
Parameters considered constant are the gravitational constant 
g, sea level standard density Po? aircraft wing area 5, 
aircraft weight Woo and the altitude of the tropopause under 


standard atmospheric conditions Hy These constants are 





32.1725 diee 


ac 
0. 7 0.002378 slugs/ft.? 

S = 400 ft.* (6.4) 
W_ = 39,000 1bs. 
H. = 36,089 ft. 


The remaining parameters in equations (6.1) thru 
(6.3) vary with flight and atmospheric conditions. These 
variations are represented in either tabular form or by 
empirical relations. These tabular and empirical relations 
are critical to the problem as they represent the mathemat- 
ical equivalents of the factors a through e listed in the 
introduction to this section. | 


It is convenient to define the constant 


(65,5) 





with the definition (6.5) incorporated the state equations 


gre 


M = gin cosa — È siny _ cagM" C, C6 6) 

y= gU ging + caoMC, ~ e (6.7) 

h= 2 siny . (6.8) 
L 
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2. Empirical Relations 


Empirical relations are used for air density ratio o, 
maximum Mach number U and the speed of sound a as functions 
of normalized altitude h. Air density ratio and normalized 


altitude are defined by 


II > 


g = L (6.9) 
Po 
and 
A 
H 
h = = (6.10) 
Hr. 


These empirical relations which are discussed in Appendix D 


are repeated here for convenience. They are 


-c,h -cah 

E c b. c4he ° Corals 

My = 2.1 - 1.1e"?:"h (6,12) 

a=a (l-c.h) , h<1 | (6.13) 

= D ccu no lc. (6.14) 

where 

e, 7 1.54100 (6.15) 

Cy 7 1.80445 (6.16) 

C2 = 0.4130 (G ET) 

Cy = 0.1331 ° (6.18) 
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3. Tabular Functions 

In situations where parameter variations cannot be 
adequately represented by empirical formulas, a table of 
values is used. Tables are used for lift and drag coeffi- 
cients as functions of Mach number and angle of attack for 
a typical supersonic fighter. eek from these tables 
are presented in Appendix B. 

The thrust Th appearing in equations (6.6) and (6.7) 
is normalized maximum thrust as it is assumed that since 
the aircraft must climb to altitude in minimum time, its 
power plant will always be operated at maximum thrust. 


Maximum thrust is normalized with respect to sea level 


static maximum thrust Thy and is given by 


O 
Thy 
JURE o l (6.19) 
M 
, O 
where 
Th, = 34,000 lbs. (6.20) 


Maximum thrust is given as a tabular function of Mach number 
and altitude for the fighter under consideration. 
4, Performance Measure 


The performance measure for this problem is 


-f de 4 (6.21) 
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5. Inequality Constraints 
The following state and control inequality constraints 


Ere imposed: 
(6.22) 
-6° =a <a<a, = 200 , (6.23) 


The maximum Mach number My constraint represents the 


"placard" limit. My 1s a function of altitude and is given 


by an empirical relation as discussed in Section VI.2. 


The angle of attack constraint a,, is set at an 


M 
angle of attack slightly above that for maximum lift coeffi- 


cient. The minimum angle of attack a, is set slightly below 


L 
that for minimum lift coefficient thus simulating aerodynamic 
stall. 

6. End Conditions 


The initial conditions are 


M(0) = Dus Ono (6.24) 
y(0) = 0 (6.25) 
h(0) = ao (6.26) 
The final condition is 
h(T) = h, = — ft. (6.27) 
L 
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B. THE EPSILON METHOD FORMULATION 
l. The Augmented Performance Measure 
Using the penalty functional described in Section 111 


for inequality constraints, the augmented performance measure 


is 
ill 
J =f at 
a Jo 
mn d 2 
+ l E E gTh cosa + & siny + caoM°c | at 
€ Jg a a D 
1 Tr. eTh sina geosy Ç 
L a= MES S 
er |» a M cacMC, i aM | at 
| (628) 
Manz 2 
+ if Ë _ = siny| dt 
0 L 
T 2K T 2K 
+r lar - 1 Pat + rf [A] P a 
0 M 0 M 


where e and r are weighting factors, 6 is a constant used to 
make minor adjustments to the admissible regions, and d Is 
the midpoint of the admissible angle of attack region; that 


is 


l + 
å = "I (6.29) 
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The required elements of w are 

















ES vU g 2 At Vs 
Wy -|M, - a A, Fe St ca o, op | | " (BO) 
| kesal, e uuu K 
a en s At s 
Wy yk = [w a M - cag kk t IE , (6.31) 
se S ,K 
a. M o 
Da thy, 2 q siny | [=] Be K (6.32) 
2M, K, T e 
» - EZ) (CE Kee isa aera kK (6.33) 
3K+k Lame J 5 x“ 5 5 
k 
0%) d K. L 
Eu" E (rat)? , pou o ER (6.31) 
Woya, = [(K-1)4t]* (6.35) 


2. Functional Expansions, Unknowns, and Partial 
Derivatives 


The state and control expansions are of the same 
form as the problem in Section V and are not shown. The 
elements of the c vector are 


T 
e = (84,85, «« (,8ys b4 Dog eee Duy Cy alan eee a Crys 


(6.36) 
dy sdas. ee , dy s Mg Yy 04 504, At) 
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where the a, 5 represent the Mach number expansion coeffi- 
cients, the b. 5 represent the vertical flight path angle 
coefficients, the Cm 5 represent the altitude coefficients, 
and the d, 's represent the angle of attack coefficients. 

The first and second derivatives of the empirical 
relations (6.11) thru (6.1!) are required. These expressions 


Bre 


=C.h -Ch : 
> | 
2 -c,h -C4h 
A p e 5 e 
Ta OT [(ez-0703h)c, + Cc, Je (6.38) 
aM 
= = 2. 61e 7? h (6.39) 
2 
d 
N 2.6.3366"? 4h (6.10) 
dh 
da . ’ 6.41) 
dha em oT PA s; 
= 0 : he (6.42) 
2 
< = 0 2 (6.43) 
dh 


The first and second partials of the tabular functions 
for lift coefficient, drag coefficient, and maximum thrust 


with respect to their independent variables and the elements 


100 


of w given by equations (6.30) thru (6.35) with respect to 
c are obtained in the same manner as in the problem in 


Section V. 


C. RESULTS 

Two problems were solved. In problem A the aircraft was 
to elimb from sea level to 60,000 feet in minimum time. 
Problem B encompassed a series of problems. The results of 
problem A were used as a first guess for the solution of a 
minimum-time-to-climb profile from sea level to 61,000 feet, 
which in turn was used as a first guess for a climb to 
62,000 feet, etc. In this manner optimal control and state 
trajectories were obtained for minimum-time-to-climb profiles 
from sea level to altitudes from 60,000 to 70,000 feet in 
thousand-foot increments. In both problems 8 coefficients 
for each expansion and 41 time points (40 time intervals) 
were used. 


1. Problem A - A Climb from 0 to 60,000 Feet 


The initial guess for the c vector was: 


- all expansion coefficients = 0 


m = jo 
Y(T) = 459 
(0) = 1° (6.448) 
a(T) = 1° 
At = 120/40 sec. (T = 2 mín.) 


Four sub-problems were required to solve the problem. 


Tables 7 and 8 summarize the performance of the algorithm. 


TOT 





Ur Om $'33" 
! 


OST x or FMMMMMMFF 3' ag" 
0.67 x 107? MPFFFFFF 3! 38" 
0.5 x Wome 2159" 





* number of iterations required 
** M — MNR method 
F — FNR method 


Table 7 


Weighting factors, optimization strategy, and C.P.U. 
time for o performance problem A. 


0.4496 x 10° 


0.2154 x 107° 


-2 


0.1544 x 10 0.1904 


-12 


ONGC TT 10 > 0.8838 x 10 





Table 8 


Components of the minimum augmented performance measure 
for climb performance problem A. 
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Neither the angle of attack or maximum Mach number 
constraints are active in this problem, The largest Mach 
number attained in the climb is 1.53 at an altitude of 
23,000 feet. Consulting Appendix D, the maximum Mach number 
at this altitude is 1.88. However, since both constrained 
parameters approach théir boundaries, a large value of K, 
is required to obtain small Jo contributions to the aug- 
mented performance measure. 

Figure 18 is a plot of the augmented performance 
measure vs, iteration number, 


: SUB-FROBLEM > 3 A 


UJ 
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iteration nunber 
Figure 18 


Augmented performance measure vs. iteration 
number for climb performance problem A. 


103 





Iterations performed by the FNR method are indicated by a 
solid line. Iterations performed by the MNR method are 
dicated by a broken line. As observed in the previous 
problem, the FNR method results in excellent terminal 
convergence. The MNR method performs well when the c 
vector ls far from optimum as indicated by relatively large 
augmented performance measures. At the commencement of 
sub-problem 1 the MNR method fails presumably because the 
initial guess for the c Vectorris too Eclose CO optimum, 

The FNR method does not reduce the augmented performance 
measure but manages to salvage the first iteration. On 
iteration number 2 the opposite occurs, The c vector is 
too far from optimum for the FNR method to work. The MNR 
method comes to the rescue. The same thing occurs at the 
beginning of sub-problem 2. In these two sub-problems the 
use of both methods in combination has allowed the algorithm 
to proceed where the exclusive use of either method by 
itself produces a divergent condition from which Phe 
algorithm cannot recover. 

Figure 19 is a plot of the performance (final time) 
vs. iteration number. 

Figures 20, 21, and 22 are plots of the states vs. 
time. In each plot two curves are shown: one is the 
expansion for the state as computed at the end of the last 
sub-problem; the second is the state trajectory obtained by 
integrating the state equations with the optimal control 


expansion, 
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Figure 19 


Performance measure vs. iteration number 
for climb performance problem A. 
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Figure 20 


Mach number vs. time for 
climb performance problem A. 
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Figure 21 


Altitude vs. time for 
climb performance problem A. 


^ integration 


expansio 


3 
Time (min.) 





Figure 22 


Vertical flight path angle vs. 
time for climb performance problem A. 
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Close observation of Figures 20, 21, and 22 reveals a 
trajectory very similar to that described in the introduction 
to this problem. The trajectory begins with a sub-sonic 
climb to an altitude of 33,000 feet. The climb angle during 
this portion of the climb reaches a maximum of 27 degrees. 
At this point an acceleration is performed to a Mach number 
of 1.53 with the aircraft ín a slight descent. A "zoom" 
climb is then performed to the desired altitude of 60,000 
feet. 

Figure 23 is a plot of the angle of attack expansion 


computed at the end of the last sub-problem. 


Angle of attack (degrees) 


3 
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Figure 23 


Angle of attack vs. time 
for climb performance problem A. 
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Initially, the angle of attack is decreasing corresponding 
to the initial acceleration of the aircraft from a starting 
Mach number of 0.6, As the eure: to climb 
attitude, the angle of attack increases. As the aircraft 
levels off for the supersonic acceleration, the angle of 
attack decreases correspondingly. The angle of attack 
begins to increase again as the "zoom" climb attitude is 
established. A further increase is evident as the aircraft 
slows down in the climb until as the final altitude is 
approached, the aircraft is near stall. The climb was 
completed in Y minutes and 18:seconds. l 
| Figure 24 is a plot of altitude vs. range where both 


quantities are obtained by integration. The scales are not 


the same. 
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Figure 2! 


Altitude vs. range for 
climb performance problem A. 
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The problem was also solved with the maximum Mach 
number restricted to 1.0 throughout the flight regime to 
obtain a comparison of the "zoom" climb technique to a 
totally subsonic climb to 60,000 feet. The aircraft was not 
able to complete the climb. The altitude of 60,000 feet is 
apparently above sap ss ass ceiling of the model, After 
H minutes and 18 seconds, which was the time required to 
complete the climb by the "zoom" technique, the aircraft 
was passing 43,000 feet and climbing very slowly. 


2. Problem B - Optimum Climbs to Altitudes 
from 60,000 to 70,000 Feet 


Table 9 depicts the results for each sub-problem as 
minimum time-to-climb profiles are generated for final 
altitudes of 60,000 to 70,000 feet in thousand-foot increments. 
For each sub-problem the results of the previous sub-problem 
were used as a first guess for the new trajectory. The 
stable behavior of the FNR method near the minimum is respon- 
sible for the ability of the algorithm to generate neighboring 
optimal trajectories. Thus, in effect, ten problems were 
solved with minimum effort by taking advantage of the results 
of each problem in turn. If, however, the change in end 
conditions is too large, the MNR method may be required to 
start the new sub-problem. This is the case in sub-problems 
11 thru 14. | 

Figure 25 is a plot of altitude vs. range where 
both quantities are obtained by numerical integration showing 


the climb trajectories for several of the final altitudes. 
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Table 9 


time to 


climb 
(9%) 
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Results of climb performance problem B. 
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Figure 25 


Altitude vs. range for 
climb performance problem B. 


The climb profiles shown in Figure 25 reveal the following 
characteristics: 

a. The sub-sonic climb profiles are identical for 
all final altitudes for the initial portion of the climb; 

b. as the final altitude increases, the altitude at 
which the aircraft levels off for the supersonic acceleration 
increases by approximately 15 to 18 percent of the final 


altitude; 


13182 





C. the aircraft performcs a diving maneuver to 
transit the transonic region with the maximum dive angle 
(13°) reached in the climb to 60,000 feet; 

d. as the final altitude increases, the aircraft 
performs the supersonic acceleration at higher altitudes 
with less altitude loss in acceleration. The results are 
not in agreement with standard practice in which the 
accelerations are generally performed in level flight at 
the tropopause. The results are in agreement with the 
results of Bryson [Ref. 2] which clearly show the dive 
associated with transiting the transonic region. Bryson's 


results were obtained by the method of steepest descent. 
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VII. AN AIR COMBAT MANEUVERING PROBLEM 


In this section the turning performance of a supersonic 
fighter is considered. First, the basic aircraft limitations 
on maneuvering flight are reviewed. Second, turning perfor- 
mance in the horizontal plane is reviewed from a theoretical 
point of view. The "corner" velocity concept familiar to 
fighter pilots is presented. Third, ane performance in 
three-dimensions is discussed. Finally, a three-dimensional 
problem is solved in which a supersonic fighter is required 
to execute a 180° course reversal in minimum time with the 


initial and final altitudes specified. 


A. THEORETICAL TURNING PERFORMANCE 
l. Aircraft Performance and Maneuvering Limitations 

A tactical fighter must be highly maneuverable. An 
important consideration in maneuverability is the ability 
of the fighter to turn. Two basic performance criteria in 
turning performance are: 

a. Ne of turn, and 

Bem Trate of turn. 
In air combat situations it is often desirable to perform 
a turn so that the aircraft's radius of turn (curvature) is 
minimized or the aircraft's rate of turn is maximized. The 
ability of a fighter to minimize radius of turn or maximize 


Beute of turn is limited by: 
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a. maximum thrust, 
bez aerodynamic stall, or 
c. maximum allowable load factor. 
2. Turns in the Horizontal Plane 
If an aircraft is restricted to move in a horizontal 
plane only, turning performance is easily analyzed. Using 
the assumptions given in Appendix A and the added assumption 


that \ 
T sin a << L, a 
ENuuations for lift coefficient, radius of turn, rate of 


turn, and the thrust required to maintain level flight at 


constant velocity are easily derived and well known. They 


are 
2W,n 
C y (T2 
L pSv 
y? 
R = nt on er) 
g(n^-1)* 
2 L 
. xm -1)* 
Y = en 9 (7.4) 
and 
To. Ey cp. MIT 
2 COS Q 


+15 | 








where 


ENTUS ta 

= angle of attack, 
REM, 
Ewlift.coeffickent, 
- aircraft weight, 


e $ 
N 


= load factor 

= alr density, 

wing area, 

= aircraft velocity, 

= gravitational constant, 
= drag coefficient, 


J 


= radius of turn, and 
- horizontal flight path angle (heading angle). 


ul Tl c E A 1 1 
li 


Figure 26 is a plot of lines of constant thrust, constant 
radius of turn, constant rate of turn, and maximum lift 


coefficient on velocity-1oad factor (V-n) diagrams. 
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Figure 26 


- Velocity-Load Factor Diagrams 
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In Figure 26(c) the constant thrust lines indicate the 

thrust required to maintain a steady state turn at a specific 
load factor and velocity. The corner velocity Ne is defined 

as that velocity at which an aircraft is capable of operating 


at maximum lift coefficient Cr, and maximum structural load 
M 


factor n, at the same time. This is the velocity which 


M 
produces minimum radius of turn and maximum rate of turn as 


can be seen from Figures 26(a) and 26(b). The corner 
velocity can only be maintained in the steady state if the 
aircraft has enough thrust available to allow the maximum 


thrust curve to pass above the corner created by the 


pm M 
not available to allow the corner velocity to be maintained 


C boundary intersection. If sufficient thrust is 
in the steady state, which is typically the case, the air- 
craft must either degrade its turning performance by moving 
off the boundary intersection until the maximum thrust curve 
is encountered, or sacrifice altitude. In this case the 
puberty for maximum rate of turn is larger than the velocity 
for minimum radius of. turn. | 

As can be seen from the previous discussion, 
optimization techniques are not required to analyze turns in the 
horizontal plane. This arena, however, is an excellent 
place to test an optimization technique which is being 
considered for use in solving more complicated problens 
involving three-dimensional maneuvering. With this in mind, 
two problems involving turns in the horizontal plane were 


solved by the epsilon method and the answers compared to 
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theoretical results. In one problem an aircraft was 
required to perform a horizontal turn with minimum radius. 
In a second problem the aircraft was required to turn 
through a given heading change in minimum time which is 
equivalent to maximizing rate of turn. The aircraft was 
given a large thrust capability so that the turns were not 
thrust limited. The mathematical model used in the problems 
is given in Appendix A. The model is an accurate three- 
degree-of-freedom model of the aircraft's motion with 
maneuvering limitations included. From the previous 
discussion the aircraft should have performed the turn 


in both cases at the corner velocity where 


muje 


M š M 


Ve ~ | SC =e imos - (7.6) 
D t 
In this equation a is the angle of attack for maximum lift 


Eerficrent and C is the drag coefficient" for maximum rift 


D 
coefficient. The a method solved both problems 
successfully. in each case the optimum trajectory involved 
wan at the corner velocity. Thus, the ability of the 
second-order epsilon method to handle problems of this type 
was demonstrated. 
3. Turns in Three Dimensions 
The analysis of turning performance in three 


dimensions is quite complicated. In this regime modern 


optimization techniques are the only method of solving 
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meaningful problems. Even optimization techniques are apt 
to have a difficult time with three-dimensional maneuvering 
problems using realistic models of the aireraft's motion 
because of large computer time and storage requirements, In 
this section the epsilon method is used to solve an impor- 
tant three-dimensional maneuvering problem often encountered 
in air combat. 

In many air combat situations a fighter pilot is 
faced with the requirement to reverse his course as rapidly 
as possible. Generally, the pilot has in mind a specific 
altitude at which he would like to complete the maneuver 
which is dictated by his desire to track an enemy aircraft 
or perform some attacking maneuver. With this in mind, the 
problem posed for solution by the epsilon method involves a 
supersonic fighter which is required to execute a 180° 
reversal in minimum time. The aircraft must begin the 
maneuver in level flight and recover in level flight at the 
entry altitude. The accepted maneuvers used to accomplish 
this task developed over years of combat experience are the 
high-speed yo-yo and the low-speed yo-yo maneuvers. If the 
aircraft begins a reversal at a flight speed higher than 
its corner velocity a high-speed yo-yo is called for and 
vice-versa. A high speed yo-yo consists of a climbing turn 
followed by a descending turn. A low speed yo-yo consists 
of a descending turn followed by a climbing turn. If the 
aircraft begins a reversal at its corner velocity, a level 


turn is called for. The purpose of applying the epsilon 
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method to,this problem is to either confirm or challenge 
the effectiveness of these experimentally developed 
maneuvers by the use of an optimization technique. The 
assumptions applied to the problem, the coordinate system 
and nomenclature used, and the derivation of the equations 


of motion are presented in Appendix A. 


B. PROBLEM FORMULATION 
1. State Equations 


The state equations derived in Appendix A are 


gTh 


M = ma Thcosa - = siny - gps MC, (7.7) 
y - & gel (7.8) 
_— a (7.9) 
h m Siny : ne) 


The states are Mach number M, horizontal flight path angle 
V, vertical flight path angle y, and normalized altitude h. 
The controls are bank angle $6, normalized thrust Th, angle 
of attack a, and load factor n. 

In addition to the four state equations, an equality 


constraint which must be satisfied is 





Th 2 
e M pša 2 _ 
0 = Wo Thsina t 2W. M Cr n ' (7.11) 
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This equation is a result of the definition of load factor 
n and is derived in Appendix A. It is possible to use 
equation (7.11) to eliminate the load factor control from 
the state equations, but this is not desirable for two 
reasons: first, the resulting state equations would be 
further complicated thus increasing the analytic workload 
required to compute first and second partial derivatives; 
second, the incorporation of the important load factor 
inequality constraint would be unnecessarily complicated. 
Parameters considered constant for this problem are 
the gravitational constant g, aircraft maximum thrust Thy» 
aircraft weight Wo? speed of sound a, base altitude Hy s air 


density p, and aircraft wing area S. These constants are 


g = 32,1725 Be 

Th, = 21,000 lbs. 

W = 39,000 lbs. 

a AD f£t./Ssec. (7.12) 
H. = 10,000 ft. 

p = 0.001756 slugs/ft. = 

S = 400 ft.? 


It is convenient to define the constants 








A | 
c, = Ë | (7.13) 
A gTh 
C57 y = (od 
e 
A 
e, = ppsa | (7.15) 


121 ' 








c m (7.10) 
lj Wo 
A 2 
E = gue . (7.17) 
e 


With these simplifying definitions incorporated the state 


equations are 


2 


M = c,Theosa - c,siny - cM Ch (TS) 
Y= c m (7.19) 
Y = c, Bone - c (7.20) 
h = T siny (7.21) 


and the additional equality constraint is 


on 2 
0 = c,Thsina + c M Cp no. (7.22) 
2. Tabular Functions 
Tables are used for lift and drag coefficient as 
functions of Mach number and angle of attack for a typical 


supersonic fighter. Excerpts from these tables are provided 


in Appendix B. 
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3, Performance Measure 


The performance measure for this problem is 


T 
J - f dt. (7.23) 


O 


4, Inequality Constraints 


The controls must satisfy 


0< Th< 1 (7.24) 
O<aca, = 24° (7.25) 
O<n En 6.5 g's (7.26) 
DI< HEN. (7.27) 


The minimum allowable normalized thrust, angle of attack, 

and load factor are approximated by zero as these constraints 
are not anticipated to be active. A zero value of the 

lower bound simplifies the associated penalty term in the 
augmented per formance measure. The maximum angle of attack 
is fixed at a value slightly higher than the angle of attack 
for maximum lift coefficient as given in the tabular data 
thus simulating aerodynamic stall. The structural load 
factor upper bound is 6.5 g's, a standard value from fighter 


aircraft operational limitations. The bank angle constraints 
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were required to keep the algorithm from generating bank 
angles greater than 1809, 
5. End Conditions 


The initial conditions are 


M(0) = 0,9 (7.28) 
y(0) = 0 (7.29) 
Y(0) = O (7.30) 
h(0) = h _ = 15,000 ft. (fast) 


O 


The final conditions are 


. v(T) 7» m Chae 
y(T) = 0 SS) 
h(T) = h. . (7.35) 


C. THE EPSILON METHOD FORMULATION 
1. The Augmented Performance Measure 
Using the penalty functionals described in Section 
BENI for the inequality constraints, the augmented perform- 


ance measure is 
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; 2 Ë 
E = c,Thcosa + c, siny + CoM Cy at 


T 
r de nsind]2 
E D È ey Mesas | dt 





T re 
dl ES ncosó eos |? 
CR pos egets, mrs 
Ir. 
+ Lf E - a siny]? dt 
E Jo EL 
(7.35) 


T 
1 š: e e 2 
+ D | c, Thsina + cgM Cr, - n| at 
T 


where € and r are weighting factors, and ô is a constant 
used to make minor adjustments to the admissible regions. 


The required elements of w are 


W. = |^. - c,Th, coso, * c, siny, + g Cp | =) des uus (7.36) 


: n, sinó 
1 F Ne I 
Wk+k D EST cmi a (7.37) 
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. cos COSY. t % 

Worn = I a o En fa » dos] NT TIME (7.38) 

E «e At 1% 
WaKtk ^ [h ü gsm ET » k=1,2,...,K (7.39) 
a At 1% " 

Wik+k = [° m sina, + eM Di n. | Ë — k71,2,. Jer E (7.840) 
ZW 1K > 

PEKIK = E = 1] P(pAt)* 5 Bee (7.41) 
20, K L 

Wer 7 E = 1 P(rAt)* , k-21,2,...,K (7.42) 
on K 
26, K a 

VOK-k ^. pes ú | ME A VK (7.44) 

Woga) = [(K-1)4t]“ | (7.45) 


2. Functional Expansions, Unknowns, and 
Partial Derivatives 


The state and control expansions are of the same 
form as in previous problems and are not shown. The elements 
of the c vector are 


T 
E (34,85, + san ,b4 Da see ,b4,,04505, «5 0y5 045055 «9 dy, 


lh. sh At) (7.846) 


Mor dy bp 1° rn e 
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where the an S represent the Mach number expansion coeffi- 
cients, the b. 3 represent the horizontal flight path 
coefficients, the SENS represent the vertical flight path 
angle coefficients, the da 5 represent the altitude coeffi- 
clents, the Em 5 represent the bank angle coefficients, the 
fa 8 represent the thrust coefficients, the Bn S represent 
the angle of attack coefficients, and the ha S represent 
the load factor coefficients. 

The first and second partial derivatives of the 
tabular functions for lift coefficient and drag coefficient 
with respect to their independent LDAP and equations 
(7.36) thru (7.45) with respect to the c vector are obtained 


in the same manner as in previous problems. 


D. RESULTS 

Three problems were solved. In problem A the aircraft 
must perform the 1809? reversal in minimum time starting 
from an initial Mach number of 0.9. In problem B the air- 
craft must perform the reversal starting from its corner 
Mach number which from equation (7.6) is 0.708. In problem 
C the aircraft must perform the reversal starting from an 
initial Mach number of 0.5. In all problems 8 coefficients 
for each expansion and 21 time points (20 time intervals) 
are used. 

1, Problem A 

Since the initial Mach number is above the corner 


Mach number, a high-speed yo-yo maneuver is called for by 
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accepted tactics. With this in mind an initial guess of 
the c vector was made to reflect this type of maneuver. 
Accordingly, the following coefficients were given non-zero 


initial values: 


a, = -0.309 
e, 7 0,524 
d; 27025533 
; TO 
e4 7 1.047 
e, = 0.262 
h; = 3,000 
The remaining initial values for the c vector were: 
remaining expansion coefficients = O 
M(T) = 0.9 
$(0) = 0° 
¢(T) = 0° 
Th(0) = 0.88 (30,000 lbs.) 
Th(T) = 0,88 (30,000 lbs.) (7.48) 
a(0) = 5° 
ACPI =_5° 
n(0) = 1.0 
n(T) = 1.0 
At = 12/20 sec. (T = 12 sec.) 


The problem was solved in six sub-problems. It took 
17.65 seconds to complete the turn. Figures 27 thru 30 
are plots of the control expansions as computed at the end 
of the last sub-problem. 

From Figure 28 it is seen that the maximum thrust 
constraint is active for the first 10 seconds of the turn. 
At t = 6 seconds there is a short period in which the 


thrust is slightly inadmissible, This is due to the use 
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Figure 27 


Bank angle vs. time for 
turning performance problem A, 
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Figure 28 


Thrust vs. time for 
turning performance problem A. 
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Figure 29 


Angle of attack vs. time 
for turning performance problem A. 
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Figure 30 


Load factor vs. time for 
turning performance problem A. 
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of the factor 6 = 1,03 in the inequality constraint penalty 
terms in equation (7.35) which has the effect of slightly 
increasing the size of the admissible region. This is 
desirable, however, as the epsilon method generates only an 
approximation to the optimal control from which the true 
optimal control must be deduced. It is easier to vr ze 
an optimal control expansion which is on a constraint 
boundary with the factor Š included, As shown in Figure 6 
in Section V, 6 - 1,03 1s the proper choice for a final 
K - 30. During the last portion of the turn, the aircraft 
is operated at the angle of attack for maximum lift coeffi- 
cient (approximately 22°). The bank angle and load factor 
constraints are not active during the maneuver. 

Figures 31 thru 34 are plots of the states vs. time. 
In each plot two curves are shown: one is the expansion for 
the states as computed at the end of the last sub-problem; 
the second is the state trajectory obtained by numerically 
integrating the state equations with the optimal control 
expansions. An observation of these plots reveals that a 
high-speed OSO maneuver is performed although the altitude 
excursions are not as large as this author expected, The 
optimization procedure settles on a nearly level hard turn 
at high load factors, steep bank angles, and maximum thrust 
for the majority of the turn. When the maximum thrust 
boundary is not active, the aircraft flies at the angle of 


attack for maximum lift coefficient. 
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Figure 31 


Mach number vs. time 
for turning performance problem A. 
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Figure 32 


Horizontal flight path angle vs. time 
for turning performance problem A. 


expansion 


integration 


Figure 33 


Vertical flight path angle vs. time 
for turning performance problem A. 
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Figure 3! 


Altitude vs. time for 
turning performance problem A. 


2. Problems B and C 

Figures 35, 36, and 37 are plots of cross-range vs. 
range, altitude vs. cross-range, and altitude vs. range 
obtained by integrating the equations of motion with the 
optimal control expansions found in problems A, B, and C. 
An observation of Figures 35, 36, and 37 reveals that the 
expected maneuvers are performed for each initial Mach 
number. In problem B the aircraft performs an essentially 
level turn from an initial Mach number equal to its corner 
Mach number at this altitude. In problem C the aircraft 
performs a low-speed yo-yo maneuver losing a maximum of 
800 feet after 90° of turn from an initial Mach number 
below its corner Mach number, In problems A and C, however, 
the aircraft does not go through as much of an altitude 


excursion as anticipated by the author. Since in fighter 
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Figure 35 


Cross-range vs. range for 
turning performance problems A, B, and C. 
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Figure 36 


Altitude vs, cross-range for 
turning performance problems A, B, amd C. 
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Figure 37 


Altitude vs. range for 
turning performance problem A, B, and C. 
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Bones, however, there are no rules on the amount of 
altitude which should be gained or lost in a yo-yo maneuver, 
a quantitative evaluation of the results is purely subjec- 
tive. The important result is that the optimization tech- 
nique did require the aircraft to perform the high-speed 
and low-speed yo-yo maneuvers predicted by accepted tactics. 
The accepted tactics are, therefore, qualitatively correct. 
The optimal times required to complete the turn for. 


each of the three problems are given in Table 10. 


Optimal time 
for reversal 


problem A 17.6 sec, ' 


problem B 20.8 sec, 


problem C 21.9 sec. 





Table 10 


d 





VIII. SUMMARY AND CONCLUSIONS 


A number of realistic problems in aircraft and missile 
performance optimization have been solved by the use of a 
second-order epsilon method. The mathematical models have 
portrayed aircraft and missile dynamics in an accurate 
manner with particular emphasis placed on the modeling of 
performance and maneuvering limitations. 

The state and control inequality constraints generated 
by these limitations have been handled by a new computa- 
tionally superior penalty functional. Three desirable 
theoretical properties of this penalty functional have been 
shown. 

A full Newton-Raphson method for minimizing the aug- 
mented performance measure:has been shown to be computation- 
ally feasible and superior in certain situations to the 
"modified" Newton-Raphson method proposed elsewhere. 

The following observations are significant with respect 
to the second-order epsilon method. 

a. The MNR method is relatively insensitive to the 
starting values of the unknowns C. The FNR method diverges 
for starting values of c which are far from optimum. 

b. Once c is close to optimum, the FNR method 
converges rapidly whereas typically the MNR method either ~ 
diverges, oscillates, or converges slowly at best. 

c. In relatively simple problems the MNR method is 


capable of obtaining a solution by itself. In more difficult 
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problems such as those solved in this dissertation, a 
combination of the two methods is required. Typically, the 
most effective procedure involves using the MNR method 
initially followed by the FNR method when successive itera- 
tions yield "small" improvements in the augmented perform- 
ance measure. In other rare cases where the initial guess 
for the c vector is close to optimum, the FNR method must 
be used initially. 

d. The power of the FNR method close to the minimum 
can be used to advantage to obtain with minimum effort 
optimal control and state trajectories for problems with 
neighboring end conditions by using the solution to a basic 
problem as a first guess for the new problem. 

The solutions to the problems solved have ‘a number of 
applications. In the missile intercept problem (Section V) 
minimum-time optimal trajectories may be used as a basis for 
comparison with the performance of more practical sub- 
optimal controllers such as proportional navigation for a 
Pont range air-to-air missile. In the three-dimensional 
turn-reversal problem the qualitative optimality of an 
experimentally developed air combat maneuver is shown for 
the first time. A significant lesson to be learned from 
the results of this problem is the importance of thrust in 
comparison to lift coefficient in the maneuvering capability 
of a fighter aircraft. Thus, an optimization method of the 


type used in this work applied to realistic performance 
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problems has application in the evaluation of tradeoffs in 


the design of future flight vehicles. 
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APPENDIX A 


MATHEMATICAL MODELS 


In this Appendix the mathematical models used in the 
problems are derived. In Section A.1 the basic equations 
of motion of an aircraft in three dimensions are derived 
under appropriate assumptions. This model is used in the 
problem solved in Section VII. In Section A.2 the aircraft 
is restricted to move in the horizontal plane only and the 
appropriate adjustments are made to the three-dimensional 
Medel. In Section A.3 the aireraft is restricted to move 
in the vertical plane only and the appropriate adjustments 
are made to the three-dimensional model. “This model is 
used in the problem solved in Section VI. In Section A.4 
the mathematical model for the missile intercept problem 
solved in Section V is derived. 


1. The Mathematical Model for an Aircraft Maneuvering 
in Three Dimensions 


In this section the basic three-degree-of-freedom 
equations of motion of an aircraft are derived. The 
assumptions are 

a. the earth is flat, 

D. the aircraft is a point mass, 

c. the mass of the aircraft is a constant, 

d. the aircraft is always in balanced flight, 

e. the aircraft rolls about its velocity vector, 


and f. acceleration due to gravity is a constant. 
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The coordinate system and notation are presented below, 





Figure 39 


Aircraft Coordinate System 


Three axis systems. are drawn in Figure 39. They are: 
a, (X AO a fixed inertial axis system; 


b. (GX! 2: 7 a monesobatine aye sea sit emgf 1xed 
to the center of mass of the 
aircraft; 


Gan UTD) a rotating axis system fixed to 
the center of mass of the aircraft; 
the x axis is oriented in the 
direction of the aircraft's velocity 
vector; the y axis points out the 
right wing. 
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The attitude of the aircraft is given by four angles as 


follows: 


a. a rotation y in the X'Y' (horizontal) plane; 


b. a rotation y in the xZ' (vertical) plane; 


C. a rotation $ in the yz plane; 


d. a rotation a in the xz plane. 


The three angles Y, Y, and $ are the Euler angles (Ref. 19). 


The angle a is the angle of attack of the aircraft using 


the thrust line as a reference. The remaining notation is 


a. forces; 


L 
D 
T 
Th 


W 
e 


l H H H I 


IA 

drag, 

thrust, 

normalized thrust, 
gross weight, 


b. angles; 


a = 


O 
ES 
(5 
cr 
(D 
un 


Ens ov 
I H H I 


d. other 


q # < 
H gw og 


angle of attack of the thrust line 
measured in the xz plane, 


angle of attack for maximum lift coefficient, 


vertical flight path angle measured in the 
xZ' plane, 


bank angle measured in the yz plane, 


horizontal flight path angle measured in 


the X'Y' plane, 


poll rate measured in the yz plane, 


pitch rate measured in the xz plane, 


yaw rate measured in the xy plane, 
angular velocity of the xyz system with 
respect to the X'Y'Z' system, 


parameters; 
velocity, 


mass, 
gravitational constant, 
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air density, 
aircraft wing area, 
lift coefficient, 


E 


= drag coefficient, 


altitude, 
normalized altitude, 


J 


radius of turn, 
Mach number, 
Corner Mach number, 


ou it 


O 


Pp Sat Die © Moo 
ton 


speed of sound, 


corner velocity, 


O 


load factor, 
unit vector (with appropriate subscript 
indicating direction), 


DS < 


e. subscripts; 


M = maximum value, 
L = minimum value. 


The equations of motion are derived following the 
methods used in Reference 19. Starting with Newton's 


second law 


dv ôv 


where = is defined as the time derivative in the xyz system. 


The aircraft velocity and acceleration may be written 


ve (A.2) 


(A.3) 


where 
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- ôv 
ve, z ét (A.B) 


and 


Meu c WU x 
s» U 


(A.5) 


US 


The angular velocity of the xyz system with respect to the 


non-rotating frame X'Y'Z' is given by 


w= pe, + qe, L (A.6) 
At this point, relations between the angular rates p, q, 
and r and the angular rates of change of the Euler angles 
are required. These relations are purely trigonometric in 


nature and are derived in Reference 19. In matrix notation 


they are 
p 1 0 -siny $ 
q = 0 cos¢ sind cosy Y (A.T) 
P .[O  -sin$ cos cosy y 


Substituting equations (A.7) into equation (A.6), we obtain 


w= ($ E psiny)e, + (ycosọ + wsind cosy)e, 
(A.8) 


+ (bcosd cosy - ysing)e,. 
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Using equations (A.8) and (A.2) the product 


u X v = (Ucosó cosy - ysing)ve, ~ (Wsind cosy + Yeos$)ve, 
(A.9) 
is formed. Isolating thrust and weight components in the 


xyz system, 


a, = Teosa, (A.10) 
Ty = 0, (A.11) 
T, = -Tsina , (A.12) 
W = -W.siny , (A.13) 
x 
W = Wi cosy sind , (A.14) 
e e 
y 
| W = W,cosy cosó , 64.15) 
Z 


equation (A.l) may be written in component form as 


Tcosa - D - W.siny = mv 5 (A.16) 


W„cosy sing mv (peosé cosy - ysing),(A.17) 


-mv (ysin$ cosy + ycos¢).(A.18) 


W cosy cos$ - Tsina - L 


The equations (A.16) thru (A.18) are the basic equations of 


motion. To apply optimization methods, it is desirable to 
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transform these equations into state variable format. 
First, lift and drag coefficients are defined by the 


expressions 


L - CL4pv'S , (A. 
D = C_2 vês (A 
D2 p ° ° 
Second, it is convenient to introduce the normalizing 
expressions 
vz , (A. 
T = Thyth : (m. 
Substituting equations (A.19) thru (A.22) into equations 
(A.16) thru (A.18) along with the expression 
Wem, (A. 
we obtain 
2 £ n ° 
ThyThcosa =4C pp (aM) N sn z aM | (A: 
Wes e | e 
W ¿cosysing = E Mípeoso cosy - ysing) (A. 
W COSY cos¢ - Thythsina _ «C, p (aM)*s (A. 


W P : 
= = aM(vsind cosy + Ycosb). 
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19) 


20) 


21) 


22) 


23) 


24) 


2559 


26) 





Solving equations (A.24) thru (A.26) for M, y, and Y yields 


the state variable format 


gTh 














M = M Theosa - & siny - RE Mc (A.27) 
W a 2W D 
a e 
e 
V E SIR Thsina sind , gpSa “e (A.28) 
Wa Mcosy UE COSY 
gTh 
D M Thsino cosó gp5a - £g cosy 
Y Wa Co wee 7 ar OS? s (A.29) 


In addition the position of the aircraft (center of mass 
location) may be required from some fixed reference point. 


To this end three additional state equations are 


X = aMcosy cosy , (A.30) 
Y = ~aMcosy siny , (A.31) 
H = aMsiny . (2.32) 


It is convenient, also, to define the load factor n as 





= TOTAL LIFTING FORCE 
n =. WEIGHT (A33) 
_ L t Tsina (A. 3) 
€ 
Th 2 
d M pSa 2 
= m Thsina + a. M Cr (A.35) 
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With equations (A.35) incorporated into equations (A.27) 
thru (A.29) the state equations are 


gTh 


iÉ LE _ (8032 2 
Wa Theosa Ssiny 2 M Ch (A.36) 





Zo 
l 


p= g nsind (A.37) 


a Mcosy 


£ g neos$ E E cosy (A. 38) 


co 


The mathematical model for the three-dimensional 
reversal problem solved in Section VII includes the state 
equations. (A.36), (A:37), (A.38), and. (A.32). Jin addition, 
equation (A.35) must be satisfied. This equation is written 
in the form 


Th 


> 
38 M pSa 2 E 
0 Ww Thsina + 5 M Cr n . (A.39) 


e 





The states are Mach number M, horizontal flight path angle 
YW, and vertical flight path angle y. The controls are 

bank angle d, normalized thrust Th, angle of attack a, 

and load factor n: The purpose of introducing load factor 
as an independent control through equation (A.35) vice using 
the state equations (A.27) thru (A.29) is to simplify the 
state equations and the incorporation of the structural load 


factor constraint. 
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The following inequality constraints are imposed 


en the controls: 


e 
|^ 


Th «1, (A.40) 
CA 
(A.12) 


The lift and drag coefficients are given as tabular 
functions of Mach number and angle of attack. Reynold's 
number effects are neglected. The parameters considered 
constant for the problem in Section VII are the gravitational 
constant g, maximum thrust Thy? aircraft gross weight Wo» 


the speed of sound a, air density p, and wind area S5. 


2. The Mathematical Model for an Aircraft Maneuvering 


in the Horizontal Plane 
In this Section the state equations (A.36) thru 
(A.38) are applied to an aircraft restricted to maneuver 


in a horizontal plane. The appropriate assumptions are 


Y = 0 (A.43) 

Y = O (A.44) 

H =H (A.45) 
O 
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Applying equations (A.43) thru (A.45) to equation (A.38), 


we obtain 


m al 
no^ ES ' (A.46) 
Substituting equations (A.43) and (A.46) into equations 


(A.36) and (A.37), we may write the state equations as 


` gu 

M = -— Thcosa = = MC, . (A.47) 
e e 

, = grand 

Y M : | (A.48) 


The mathematical model for the two dimensional 
minimum time and minimum radius of turn problems referred 
to in Section VII includes the state equations (A.47) and 
(A.48). In addition, equation (A.35) must be satisfied. 
Using equation (A.46), equation (A.35) is written in the 
form 


En 


O = uU Thsina cosd + 
e | 


2 
pSa 2 
2H CM cosó - 1 . (ma) 





The states are Mach number M, and horizontal flight path 
angle Y. The controls are bank angle ç$, angle of attack a, 
and thrust Th. It is possible to eliminate one control by 
substituting equation (A.49) into equation (A.48). The use 
St Equation (A.49) as an additional equality constraint is 


preferred, however, because the state equations are simpler 


D 





and the required control inequality constraints are simpler 
to incorporate. 
The following inequality constraints are imposed on 


the controls: 


GNE «185 50) 
2 -1f 1. 

0« ç$ < bu = cos ind ; CAG) 

0 < Qa al (A.52) 


The lift and drag coefficients are given as tabular 
functions of Mach number and angle of attack. The parameters 
considered constant for the problem are the same as those 
listed for the three dimensional model described in Section 
ES. 


3. The Mathematical Model for an Aircraft 
Maneuvering in the Vertical Plane 


In this section the state equations (A.27) thru (A.29) 
are applied to an aircraft restricted to maneuver in the 


vertical plane. The appropriate assumptions are 


q = 0 (A.53) 


and 


UO (A.54) 
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Substituting equations (A.53) and (A.54) into equations 


(A.28) and (A.29), we may write the state equations as 








E | ug . £p5a y? 
M Wa Theosa - 2 siny Zu MC, (A.55) 
wi sihy Thsina , £058 wo _ Ë cosy (A.56) 
Y Wa M aW. L a M , 


The mathematical model for the two dimensional climb 
performance problem solved in Section VI includes the state 
equations (A.55) and (A.56) along with state equation (A.32). 
It is convenient, however, to introduce the following 


relations into the state equations: 
g = E (257) 
O 
Th = — (A58) 


h= m (A.59) 


where 0 = density ratio, 
p standard sea level density, 
Th = normalized maximum thrust, 
H, = tropopause altitude, 


and h = normalized altitude. 
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SUIS Ito equations (AW57) thru (A.99) into equations 


(A.55) and (A.56), we obtain the revised state equations 


gp sa 








M mm cosa = sina u oM Co (A.60) 
š gp_Sa 

E gTh sina O _ gcosy 
E cv Ww 2W. omc, aM (4.61) 
nem i SIME Y ae (A.62) 


The states are Mach number M and vertical flight path 
angle y. The control is angle of attack a. 
The following inequality constraints are imposed on 


the states and controls: 


0 «M«M (4.63) 


Amin £ 0% < Uy : (A.61) 

Thrust (Th) represents normalized maximum thrust for 
the problem in Section VI. This is given as a tabular 
function of Mach number and altitude. The lift and drag 
coefficients are given as tabular functions of Mach number 
and angle of attack. 

Empirical relations are used for density ratio o, 
speed of sound a, and maximum Mach number My as functions of 


eae. These relations are given in Appendix D. 
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The parameters considered constant for the problem 
are the gravitational constant g, standard sea level density 
Po? wing area S, gross weight NR > and tropopause altitude 


Hy. 


4. The Mathematical Model for the Missile Intercept Problem 


In this section the mathematical model for the missile 
intercept problem solved in Section V is derived. An air- 
to-air missile launched from an attacking aircraft must 
intercept a constant-velocity target. The missile is 
restricted to maneuver in a plane. The orientation of this 
maneuver plane in three dimensional space is defined at 
launch as the plane containing the position of the missile 
at launch, the position of the target at launch, and the 
velocity vector of the target. 

The assumptions applied to the problem include those 
presented in Section A.l plus the following: 

a. iche initial velocity vector of the missile lies in 
the maneuver plane, 

b. the attacking aircraft is tracking the target at 
missile launch so that the missile's initial velocity points 
at the target at t = 0, 

c. the target es with constant velocity, 

d. components of out of plane forces perpendicular to 


the maneuver plane are ignored. 
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Figure 40 presents a view of the problem in the maneuver 


plane. 


TGT Y 





Figure 40 


Missile Coordinate System 
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Four axis systems are drawn in Figure HO. They are: 


(AS Y") a fixed inertial axis system in the maneuver 
. plane with the origin at the missile launch 

podmes 
p OC SY) a Newtonian reference system in the maneuver 


plane with the origin at the missile launch 
point at t = 0; after launch the origin 
remains fixed with respect to the target 
(it moves with velocity Vin with respect to 
the X'Y' system); 


Ce (xy) a non-rotating axis system fixed to the 
missile center of mass; 


EM (Xoy) a rotating axis system fixed to the missile 
center of mass; the x' axis is oriented in 
the direction of the missile's velocity 
vector. 

The systems X'Y', XY, and x"y" are oriented in the maneuver 
plane so that the axes O'X', OX, and cx" form the intersection 
of the maneuver plane and a horizontal plane. These axes are 
chosemmsowthat ithe tamget's imitšal Xe X*, ande positions 
are positive. The axes O!Y!, OY, and cy! are chosen so that 
the component ofmissile weight inthe maneuver plane is acting 
in the negative Y', Y, or y" direction. A11 angles are 
positive as they are shown in Figure 40 in the counter- 


clockwise direction. The remaining notation is: 


a. forces; 


N = normal aerodynamic force, 

A = axial aerodynamic force, 

T = thrust, 

Th = normalized thrust, 

W. = component of missile weight in the maneuver 


plane, 
b. angles; 
angle of attack 


missile flight path angle, 
target flight path angle, 


a 
0 
Y 
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c. other quantities; 


€ QO 0O = 
IS SE 


DS o oo os 


missile velocity, 

target ye LOGE, 

missile Mach number, 
target Mach number, 
missile gross weight, 
normal force coefficient, 
axial force coefficient, 


angular velocity of the x'y' system with 
respect to the x"y" system, 


missile mass, 
gravitational constant, 
missile wing area, 

air- density, 

load factor, 

Speed of sound, 


The equations of motion are derived following the methods 


used in Section A.l. Equations (A.1) thru (A.5) are 


identical. 


The angular velocity of the x'y' system with 


respect to the non-rotating frame x"y" is given by 


The product 


is formed. 


4=0€, (A265) 
6 x v = v O ey (A.66) 


summing forces in the x' and y' directions, we 


obtain from equation (A.1) 


ND 





mv , (A.67) 


Teosa - Nsina - Acosa - W sine 


Tsina + Ncosa - Asina - W. cose mv 6o . (2658) 


Axial and normal force coefficients are defined by the 


expressions 
J. 2 
A = C. 2 pov S 5 (A.69) 
and 
= 1 2 
N 7 C4 5 pvfS . (A.70) 


Substituting equations (A.21), (A.22), and (A.23) into 
equations (A.69) and (A.70) and transforming the results 


into state variable format, the state equations become 


gW 














"= gy E ae _ gosa c 

M aW Thcosa AU M C cosa SW Mé C sino AW sine , (ATI) 
e e € e 

- SINN Thsina _ gpSa gosa P N. cos® 

8 aW. M AL MC „Sina + a MC, cosa aW M (A.72) 


Two additional state equations are required to impose end 
conditions on the relative positions of the missile and 


target in the optimization procedure. They are 


nie 
II 


aMcos0 - aM cosy (A.73) 


T 


Sinn, (A.74) 


H£ o 
ll 


aMsino -— aM 
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The states are missile Mach number M and missile flight 
path angle 0. The control is missile angle of attack a. 


Normallzed thrust is given by 


Th 


N 
= 
ct 

|^ 
ct 


B (A.75) 


Th 


tl 
O 
ct 
V 
ct 


where Er represents engine burnout. The following inequality 


constraints are imposed: 


-Q 


|^ 


a < w (A.77) 


a 
8 


—I3 


| ^ 


(6M) < ny (A.78) 


Equation (A.78) represents a structural load factor limit. 
The axial and normal force coefficients are given as 

tabular functions of Mach number and angle of attack. 

Parameters considered constant for the problem are the 


gravitational constant g, maximum thrust Th the speed of 


M? 
sound a, missile weight Wo? air density p, missile wing area 


S, the Mach number of the target M and the flight path 


p? 
angle of the target Y. 

In order to properly define the problem, it is necessary 
to perform several manipulations in analytic geometry. 


First, the three dimensional positions of the missile and 


target must be specified at launch. Second, the velocity 
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vector of the target and the Mach number of the missile at 
launch must be specified. Once this is done it is necessary 
to: 

sa. identify the maneuwer plane, 

b. identify the XY coordinate system, 

c. calculate the target coordinates in that system, 

d. calculate the target flight path angle Y, the initial 
missile flight path angle 0(0), and the component of the 
missile weight acting in the maneuver plane Wo The 
optimization procedure can then be commenced. 

To accomplish these calculations an initial coordinate 
system i established in which the problem can be easily 
ee alized. The origin is situated at the missile. Tre 0X 
axis is positioned in the horizontal plane. The OY axis is 
positioned in the horizontal plane such that the target has 
nom, coordinate, The 0% axis 1s positioned in the vertical 
plane such that a target which has an altitude advantage 
over the missile has a positive Z component. This coordinate 
system is shown in figures 41 and 42. The angles Sm and By 
are defined as shown above. The following relations may be 


written: 


a = Roi + Ak (A.79) 

Mp = maicmaj mak (A.80) 
m, , 

tand, = — QUE RO 
Mo 
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"T 


measured in the 
vertical plane 
containing M+ 





Figure 41 


M TGT 





Figure 42 


Initial Missile Coordinate System 





y' f (A.82) 


With the problem defined in the coordinate system 
shown in Figures 41 and 42 it is now necessary to transfer 
the problem to the coordinate system used in the optimization 
procedure. That is, it is necessary to identify the maneuver 
plane and the XY coordinate system. To this end, a — 


normal to the maneuver plane is 


N = a x M (A.83) 


-h mori + (hm, ,-Rym, ,)J + RM K x (A.84) 


T e y 
To establish the X axis a vector is required which is in 
both the maneuver plane and a horizontal plane. Such a 


vector is 


X = N x k (A.85) 


(him, , -Rum, , )1 t hg, J (A.86) 
To establish the Y axis a vector ls required which is in 
the maneuver plane and perpendicular to the X axis. Such 


a vector is 
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Y= N x X (A.87) 
Z 2 
= -Rphpm, , i + Ry, , (hom, ,-Rom,,)j (4.88) 


=[hnm ^* (ham, ,-Rym, ,)* Jk 


y: 


The angle $ between the maneuver plane and a horizontal 


plane is required and is given by 


N + k 
cosq = |—— (A.89) 
iN] [kx] 
Ben 
= —— ° (A.90) 
[ (ng, , + (hm, ,-Ryn, ,) + (Rem) J 


The missile weight component in the maneuver plane We may 


be found from 


W = wsing . (Ano) 


This is shown graphically in Figure 43. 


maneuver plane 





horizontal plane 


Figure 43 
Missile Weight Component 
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The initial target coordinates are 


Xm(0) = |PROJ,a| 
N Ry (hum. ,-Ryn, ,) 
Mu ccc mE. compres (A.93) 
[ (hum. , -Rym, , ) + (hym.,) ] 
Y (0) = +PROJ ya 
T . (A.95) 


2 2 3 2 2 
uf “Rn hgm,, Ap m, -hg (hum. , -Ran, , ) 


dye a 2-915 
[run Rm, Qm sf A n aR) I. 


The sign of Yn(0) is resolved by: 


if Am > 0, Ir is positive; 


jen ho EDT Yo is negative. 


The initial missile flight path angle 6(0) is 


_ Yn(0) | 
tane(0) X (0) (A.96) 
M 


by assumption b at the beginning of this section. Before 
proceeding it is necessary to insure that the X and Y 


vectors given in equations (A.86) and (A.88) have the 
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correct sense, This may be checked by observing the sign 


of PROJ,,a which should be positive and the sign of PROJya 


X 
which should be positive if hm > 0 or negative if h 


"g 


T < ¡Qu 
After the senses of these vectors have been checked and 
altered as required, the target flight path angle Y may be 


calculated by 


x 
osy = — ———— , (A297) 
IMs XI 
The possible range of v is 
eye. (A.98) 
COS y 1s positive, then 
E PES EN. (A.99) 
MEC osy is negative, then 
s Y < T or We eS Ye See E . (A.100) 


To find which inequality applies in (A.100) the quantity 





Mn! : 
k = ——— (A.101) 
[Ma | TY | 
ioeetermed., If k is positive, then 
0 < Y < T ; (A.102) 
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If k is negative, then 
-T < Y < 0 . (A.103) 


This logic completes the set up of the problem in the 


maneuver plane, 
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APPENDIX B 


TABULAR FUNCTIONS 


In this Appendix the tabular functions used in the 


problems are presented. 


l. Three Dimensional Plots 
Three dimensional plots of all tabular functions are 


presented here. Figure 44 gives the lift coefficient C, as 


L 
a function of Mach number M and angle of attack a for the 
the supersonic fighter aircraft used in the aircraft 


problems. Figure 45 gives the drag coefficient C. as a 


D 
function of Mach number M and angle of attack a for the same 
fighter. Figure 46 gives normalized maximum thrust Th as a 

function of Mach number M and altitude H for the supersonic 

aircraft performing the minimum-time climb in the problem in 
section VI. PRE L7 gives the normal force coefficient Cn 
as a function of Mach number M and angle of attack a for the 


air-to-air missile used in the problem in Section V. Figure 


48 gives the axial force coefficient Ca for this missile. 


2. Tables 
Following each plot a condensed version of the table used 


in the computation is presented. 
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TABLE 11 (continued) 
LIFT COEFFICIENT FOR A SUPERSONIC AIRCRAFT 
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TABLE 11 (continued) 
LIFT COEFFICIENT FOR A SUPERSONIC AIRCRAFT 


MACH NUMBER 
QF ATTACK 
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TABLE 12 
DRAG COEFFICIENT FOR A SUPERSONIC AIRCRAFT 
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TABLE 12 (Continued) 
DRAG CDEFFICIENT FOR A SUPERSONIC AIRCRAFT 
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TABLE 12 (Continued) 
DRAG COEFFICIENT FOR A SUPERSONIC AIRCRAFT 
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TABLE 13 
MAXIMUM THRUST FOR A SUPERSONIC AIRCRAFT 
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TABLE 13 (continued) 


FOR A SUPERSONIC AIRCRAFT 
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RAMETER = MACH NUMBER 
PARAMETER = ALTITUDE 
1.200 1.300 


s PA 


HORI 
VERTICAL 


1.500 1.600 1.700 


1.400 


1.000 1.100 


0.900 


Q UYA OOO + O O (YAP O -OOO T à «r0 9 — CN rS Or ed 
Q C9 O OY += rJ Q QO OO 4 O QOO rF- 0 OY 00 ++ O Cm 
|D — cO O .O Of 0s 700100 030 iu NO OO F- Ou «f 0 mO 
NADOMm OLA + f$ CQrU €N ———=— OC OOO OQOCO 
ve @ @ @ @ 0 o ù @ o @ O ‘o @ o 8 0 eo 0 @ 9% 9 @ @ o @  @ ‘GH 


-AHOOOOOOOOOOOO2OO00000009 


O A $MN YD INM DML DUI CO Ae Q un (N —+ O 
O Y DMA «0 003 ++ (ACO NWONTMTTIODNNNOMm 
FA TMHMMOMAMm AHAOIO AID INO OF ansrmainetO 
NOOWMm iM I: N= ——-=—— OOO O O O OO O 
e» © @ @ 9? ọọ o 9? 0 @  @ @0 @ 0 0 9? 9? 0 0% 0 0 9? 9 0 @ 


"“00002000000000000000000 


OP-NDOMMIOTPrMIWAITO IT TONSONm 
COO* OO 00 $ 70-00 0U «f OO eQOnr- com o0vo:o 
ON HO O CN uy mU exp p. «reO]GoO0 WQUu V ANA O 
"0 DDR DUI FIMONN A —à—— O O C O OO O 


@ è o ọọ ¢ où 9 90 € 0 9 *. 9 9 * 0 Q 0. @ @ ce @ 9? e 


ede40O0O0OO0OOOOOO'DOOOOOOOOOOoooo 


ODAD (0 e c0Y0 OU IL NO U^ CJ eU CNUA CN SO e UA 
SO FO m0 ft- 40 00 20 7^0 00 7O -F. O COO enu Oreaony9 
ADP OM VO MOA PO GNALA GNAT AN HMO 
AN Or .QUALO XT CON CON A 0 10 COCO OOO O O OCO O 


ciOOOqO]QOOOOOOOOOOOOOOOOO02D00o2 


OF YO MD DN ya QO0 ent -daO-000N 
ODDO OM DO ON AMO NA DAM OM AO 
+ ry FODATDOARDOMADM OMS MN Nr O 
QOO O Ura + + QN N—4 ——+—Ç O CO O CO O O CO O O 


-AOOOOOOOOOODODOOSOODIOOO90O00900O09 


ONNOITNHATOTMIMIANNMNOOMNNDT ANNO 
O Y DA NADANDO Y PIN Dr NIT DO 
MAOAAIT OO MM NDNANO DEM OF FONAAO 
ODO D A FIONA HO OODOVDOODOOD 
@ e 9 09 e. 9 0 e@ 9? e @ 9? 0 @ e 9 9€, 0 0 @ 0 e o @ €? e$ 4 


Qe f-i-00 C0 NM ADD UAM NA FO MONA A 
DOUAMONOND TAN DONTNNDMONLOON 
Qo P- f- C0 O^ eA Xp 00 C440 000 «fF 69) OO 00 GO Un qr i+ NANN DO 
DO OMAP Q eoeJoN---AROOOOOÀOOOOOO 
t e e «9€ e * 9 9 9$ o 9 9 090 * 9*9 9? 0 0 9 9 9? @ 9? 0 e. 4 


OOOOOOOOOOOOoOoOooooooooooo 


QUO ruo 20 c0 (0 HP CUN PU CN 6900 O co oa YO «f un 
- NO O O (Gun «X cr 0YU0-F OD 00-00 C0 C00 Ou 
STOO DO YM AN AtEm E ADO DANNA AO 
Q*COP- «OUS YmÓóooN-—-—OOOoOcOoOoOoOoQoo 
* e 9 @ 8 9 9? 9 * 9 @  @ oe 9? 9? e 06 0@ @ 9$ 9 ? 


OOOOOOOOOOOOOOOOOooooooooo 


QOO OT 0400 0 --T O UA cOUYONONS T uOJ4IAMTfO.r 
+O! NS O OO —— O. OO CO UYA NA SO FFAA DONAN 
CJ CON CJ 000 CO YN DO +O XQ C eU CO Qu T eO C1 09 A O 
Q^" CO P «DU «fT On040J 7-0 OC»0O0OQO20O0QO0 
€ e 9 e @ @ 0 9*9 9 e 9* c. 0 9 @ 9? 9 0 9? 040 9 9$ @ @ @ 


QOOOQOQOOOOOOOQOQO*cOOOOOOoOOOo 


QOO OOQOOOOOOOOOOoOOoooOooooodoco 
ee 9 e 9* 9 e 9 9* 9 9? $9 9 0 ce 9 0. 90 9? o @ 0? e 9? @ 
OQ OOOOOOOOOOOOdoOOoooooooooo 
O O C O O m S lejsipleljelsiojejesibliala 
OOOOOCOCI0)0oQCOQUOOoooccoooooo 
NOWONONONINONONOWNDNONOWO 
e AN JC OS PUVLAO OP DOADARPOOAAN 

st HAs 


279 








TABLE 13 (continued) 
MAXIMUM THRUST FOR A SUPERSCNIC AIRCRAFT 
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TABLE 14 
COEFFICIENT FOR AN AIR-TO-AIR MISSILE 
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TABLE 14 (continued) 


NORMAL COEFFICIENT FOR 


AN AIR-TC-AIR MISSILE 
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TABLE 14 (continued) 
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CLE OF ATTACK 


CH NUMBER 
100.000 


ETER = AN 
ETER = MA 


M 
M 
50.000 


A 


AL PARA 
PAR 


120.000 


110.060 


70.000 80.000 


60.000 


Q (YO (YO Q^ 00 (0 O O m + O «0 400 O TAI OWANZNTOOIO AITNOWTHNODADO 
Q ec O .O O + IS 22 O MUYNUUdYOONO mo r-uYoo $ MODNDGO YA —2.n O + GO OO + O 
OO BO MLAB DOO FAP PND ADD DR FO MES FM OD OM HPP OMA OP O yg 
ANY ODBDAIOD AFP Re OPM AR OP ONTO N AOR D OF ONTONAODNODM™ O 
9€ e 9 e 9 9 ọọ 9* 9 ọọ 9" o o O o 9" 9 9 9" O o O o 9" 9 9 9? 59 c. Oo 9, 9? 09 9, 9 9 9$ 9$ 9? e? 
GO CO O OOO QY GA OOO CJ oy T T Jr xMWrsr rrr ooxmmeuwmnmameeeeec 

ten 


(MORA DONA T uYO CO 6JO CO 0J OO eo T OO Q^ UY e Jf (O cñni— — f — — O QG + — ++? OO nM 
— 5 8r J GYA O ++ (N -O O +O ni DUEVh FONDA DS SCA DOF TION > + cO N.O O 
NY MO NM O 0 T CJUO m CO SD O T QN P O t= Cs CN ° CN O r try CN Ç" r> < O DOY 
VOM 0 O (N G: CJ Y 00 P 00 UN GO F-- 4 UYCOY CJ O US CO F- .O "QUAM SMN OD U^ f OUS a CO (eO y 
e E o 9 9 ó 9? $9 € 9? o 9 9 9 96 9 0 0 0 $9 90 9 09 9 o 9 9 $9 9 9 6 9 be N 9? n n) 9 0 4 09 9 
OQ U* OQOQO Ar OI Y UND pa Pa fo Pa Pro ps Pu e Pe Pa DD OD) DD OD DD UVLA LALA LA LALA UA 

rir re Al A A A A A Al AA A A A A A AAA th dh Hh ht — — — — — — — — — — — — mM — — 


n O 9700 0* — O uA O C —:O m 9 O +f —@ .O + rn try Q Yt muy O0 OS. COU O 2 O — O O .O + -O— 
DIN Y $e ON N ONMO MANDON NMN +$ O Oun ON try uA CN O O in O rN Dio 
Ly QS — r + O F= uy t (r FDOMMH DON DMD Pe $ 2 CN cy ++ (ry CN CN A — [= MAM ODDO 
CO C^ CJt'YO $ OOLFM ONMM ADDO OWMNNQNOD OF ONS MAO OTAMAN AOO 
a. @ ç @ ©  E 1 9 9? 959» 9 9 9 9 o HH 9 9 0 0 ese 9 9" 5.» 9. 9? 5. 9 9 @ @ @ @ 9? @ 9? 0... 0 @ 
QOO ——^0OJ (f UNIO OON ONOOO OON Q^ 00 0900 0000 00 c0 00 00 P P7 P P9 P P P Pr 
teen 


— QO @ <D co (N I— 4^ 09 00 eo P OW [S CN 4 CO UY XP 09 OO h- (040 00 ^9 90 OS (4 OO «FC Q^ P CO (OU P 
— (N O iY OO CO O G O OQ rq Q OD ed O c d CJUO OO r- Og coP- OO O1unh€Og cfe CN O 
MAPPED OM AMDNNAMMODTSCINDONO OAODTPNOMALSAPOVMYUN nun 
NFOAMOLFADM OP ANH ONS NADDONS ANAODEMALTMADAOMYD 
€ 9 @ ç 9* Ó 0 9 e. 9 9 o. 9 9 9 9 e. 9 @ $9 9 $9 @ Q 0 n © @ 0? 0 @ 9* 9 @e © @ o 090 @ © @ @ 
e c0 — — (NJ (NJ Q + LY TO — 2 O O O O O Ó @ QO — O. C; Q: C" O: O° O° C*° c0 <o ñ Q G O m m 
Set — — — A A A — —— (NJ (N (N CNJ CN (N (NJ (N IN rl Al Al A Al A — — — —— —— —— — eA —— — 


OY FONNNOMO —2—4 (N ++ rN > (n Ə =— t IS + O — cJO5 rn O G*Y ++ OO -r OO CJ UY XO CN STAND 
MID ODPADMORDOFIAMA OHO ON FNANAM O O €N 92 9 OO O^ ñ — DINO TD iO 
QC UM QOSUM OO (NJ uy Q un —@ CN rN. +$ trY IS + O O (A — uy tX 9 CO CO 1 YACNIL 239 (0 0 (e 7O P COON (NJ 
DO MONO PASO OY OCNOJeQO QOO Ooh O un E Y — O —T= uy ++ n CN cN 
@ @ ç 9 9 @ 9 9 9 9 e o o 5. 9 9 0  @ 94 9 9 5 9? 9 8 9? o 9 9 o 9 b @ QQ Q @ @ @ 0... 0 
Q - —AUNOJOY en +f t I— O 2 O O O O c°2 Q, C: O°. O° G: >, Ot c5 O O O O Q cO OQ <-> Pa po Poo Pr fa Pr 
A A A — —— — (N VN (N (JN NR —— — — — — — — — — — — — — — —— —— —— — 


O00 «DO OO Uu XP e OD «T M +-O 'N O nO oJ 05 Oruun += CJ 00 e (040 ^0 4r UY CJ 09 P — (Na O O O 93 
QN < —Oury (NOOO O SN NIYA CAO a OI CY R.O T+ CG à= Oxa (rN Gn UN COS e 0NUM (0 090 00 e P o ea cy e 
LyAOO O O o (N: (y O F— Or OS OO PUO 00-9 >x O (NJ O QQ Q` — “N + (N C "O (N Ç GO Q — — N + = 
NN QOoOh- OO 79 c NOJ«O00 20 CO P O NY ALO CNT OY DMA SIMAO DA DUNAS MN NS A O 
e è 9 9 ç @ 9 ó 9* 9? ç 0 0 @ 8 9? o 0 0 Q 0 9$ 8 9 06 060 6 0 0 06 0 6 4 0 06 4 0 . 0... HF c @ 
OO OO-MoJejm «d «Or-9)00990 09005 0 ^009 00 CO P Ie f IS PP PIS DS 9.9 00000030 
ten 


LAGO O (rN (N O i— + O Gr gg AJOS SM MORAS NA) OD NO TO Am LAN FAT MO 
SQ e x c4 0 09 00 OAS MODEM NN rU ONAN Y O CJ 9 DH FAD FAMN MDN ADN N 
C 4F OO jy T (rN O O IZ< 00^. Qn 0s em UNST QUY) AUSGOS -CY:990 I< X TONO DD OO LA 
Ce 300 Uu CO O00 e ST HO p CO QN O^ OS PO P SO 4 CO NOOO UON 0 CO - OUS F 0 CJ e 8 OO P D ry T 
e 9 ọ 9* 9 o sv 9$; @ 9 8 9* @ oo? O $9 9 9 9 9* 9 Q? 9 $9 9 9 e$. 9 9 9 90 $9 c. 9? Q? 9 9 0$ ? c c 
AOD WDD © aA AN OY FUN OO OO DO OO OQ O O IO AUA UY N UYAUAUN UA: OY UA gr T Tr HT T HT HT 

A lA A A e A A A AA AAA A AA A AA A A A A A rl AA A A A A A A Ad rl A Al dl Al rl 


Q e (CONJ e xr uo P 00 0v O eC CÓ Fr USO PS O0 OY O eO 0 ST UY DD DO ANNO SN DM ODO 


QOO OO OO OQO O M —— — — — Mo o —OJJ CJ 0g 0J CJ o OJCNV CM en 00 09 00 00 0e e) 09 e r SS 


185 









YO 


SS 
LIA 
U^ < Z 
7 
=Z 

M 


if 


// 
IN / 
H w /) A 77 Z 





TABLE 15 
AXIAL COEFFICIENT FCR AN AIR-TO-AIR MISSILE 
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TABLE 15 (continued) 
AXIAL COEFFICIENT FOR AN AIR-TO-AIR MISSILE 
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TABLE 15 (continued) 


AXIAL COEFFICIENT FOR AN AIR-TO-AIR MISSILE 
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APPENDIX C 


INTERPOLATION 


In the optimization problems solved herein the aero- 
dynamic data is given in tabular form. The dependent 
variable D is given as a function of two independent 
variables M and a in all cases. Excerpts from the tables 
used in computation are presented in Appendix B. 

9D 9D ð D ð D 


For a given M and a quantities D, a) Sh > De sD» 


34D 


dMda 
Parabolic interpolation is used to obtain these quantities. 





and are required by the optimization algorithm. 
In this Appendix parabolic interpolation for two independent 


variables is derived. 


i Farabolic Interpolation in the Elane 


To apply parabolic interpolation to a tabular function 
of two independent variables the nearest point in the 
tables to the given point (M,a) must first be found. It is 
assumed that the tabular data is given at constant intervals 
AM and Aa in the independent variables. The nearest point 
given in the tables and the surrounding eight points are 
required in the interpolation and are shown in Figure 49. 
The parameters 0 and $ locate the point (M,a) from the 
nearest tabular point (M. , a De. A) is the nearest 


S 


point then 
- > < $ < Š (C.1) 
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s+1*%s+1) 


s+] s) 


5+1 19.1) 


Figure 49 


and 


GG 


I 
hole 
|^ 
O 
|^ 
hole 


The inequalities (C.1) and (C.2) hold unless the nearest 
point (M sa.) is on a border of the table. In this case 
(M a) is chosen one point in from the border. The 


parameters 8 and $ then satisfy 


- INS 9 < 1 (C. 


and 
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n2» 


3) 


ü) 





Writing Taylor series expansions including up to second 


order terms for the eight points surrounding (M. ,a..), we 









































obtain 
D(M_+AM,a_+Aa) = D(M. 14 0,44) = D(M 5a.) (CaS) 
BD 2 2 o p. - 
9D 9D (AM) 39D Gio) Sa oD o D 
G + an Hr, ro | 5 ee 
oM Š 9a a 2 amc E 2 CAP Md o E 
D(M, -AM,a -40) = D(M, 4,0, i) s D(M, 4.) (CTO) 
2 ne Qoae p 
S S 2 aM* '!s OW "S ' s 
D(M.-AM,a *Aa) » D(M. ,,a.,,) ? D(M,,o.) (0275 
2 ne Q xd 2 
9D 9D (AM) 97D (AG) Bap 2 o D 
- AM — + AG — + -= Zn De A =- AMAQ 
oM E a 2 2 amé , 2 B E MIA x 
D(M +AM,a -Aa) = D(M, 11,9, 1) P D(M, ,a. ) (C8) 
2 ne 2 ne 2 
9D 9D (AM) d=D (Aa) Ə D oD 
+ AM - Aa = + == + — - AMAa Ama 
9M S ls 2 ame s 2 Jal s mag S 
e ae 
3 9D (AM) 9 D 
= x + A —— ZZ 
D(M_+AM,a_) D(M 47>9,) D(M, ‚a, ) M SM ' - 
(C.9) 
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, = x £ aD (AM) 9°D 

Dee OS DA 7 99.) D(M 50.) - AM oo} + 2 
S oM 's 
(S10) 

D(M sa +a) = D(M ea aa SD ¿600% 2D 

s? Us a co eg g?g Q 3a 5 
S Ja S 
Comm) 
2 „2 

a x | dD (Aaj” 9D 

D(M, »a,-d0) = D(M »a,_7) D(M sa, ) - Aa Sm 4 = go 
S Ja ‘Ss 
(C .12) 


Subtracting equation (C.10) from equation (C.9), we obtain 


ES ow 
Eos c PUR TAN sy | 


or 


2 D(M. 1,0.) = D(M. 4,0.) 


SAM (0.13) 


V| 
=Ë 


Subtracting equation (C.12) from equation (C.11), we obtain 


T oD 
DM 0341) = D(M.,a. i) 2^a Ju I 
or 
3D D(M_ >. 43) E D(M.,a. ) — 
da 2^a i ° 
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Adding equation (C.10) to equation (C.9), we obtain 


2 
E 29 D 
DIM 432%) + D(M,-,,0,) * 2D(M,,0,) + (aM) > 
oM 's 
or 
3 ^D - D(MS 41 29s) E 2D(M, 20.) E D(Ms ts? (C.15) 
2 EA AA : 
9M Is (AM) 
Adding equation (C.12) to equation (C.11), we obtain 
D(M ) 4 D(M ) * 2D(M ) + (A ye ə “p 
s* 3 s+1 S sj s%g = da? a 
or 
2 2 y 
9a 's (^a) 


Subtracting equation (C.7) plus equation (C.8) from equation 


(C.5) plus equation (C.6), we obtain 


PD 


Mo-1>%5+41? 


a 3) ~ BC A 1) 
AMA |. O 
or 
^D - DM 47 20541) + D(M. 4,9. )) - D(M, ,,0,,4) - D(M f, sos 3) 
2M3 0 i AMAa 
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A Taylor series is written for the point (M,a). Using 
Sailations (0.13), C.14), (C.15), (C.16), and (C.17), we may 


reduce this series to 
š z 9 d 
D(M.*6AM,a. t$Aa) 7 D(M,a) * D(M.,o.) * 5 [D(M. ,,,o.) - D(M. 4,a.)) 


2 
$ = ° a 
+5 [D(M_ 0.47 ? A +5 [D(M 150.) 2D(M. ,a. )*D(M. .,o.)] 


(CB) 


2 
it $- [D(M, 0,41 )#D(M,,a,_, -2D(M, 501.) 


iO š E 
+ SP EDO. 9947 DO. 30 1) D(M. 1,95,442-D(M. 7 2 70) 


Rearranging, we have 


+ 


$9071) p(w ‚a, 


D(M,a) * T D(M. 1595.3) S 1? 


0 8(8+1) 
s: s: D(M. 11595 1) * 2 D(M 41 >95) 


we 
+ 


Ü (6*1) 
E D(M 1595447 NT D(M, , 0.41) (C.19) 


+ 


6(6-1) 
ge s D(M 


6 
_ T D(M. 159543) 5-1: 9g 


+ (1-0%4%) D(M_ a.) 
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Equation (C.19) is the expression used to interpolate for 


the value of D in terms of the surrounding none tabular 


points. To obtain expressions for the required partial 


derivatives, observe that 


M = M. + BAM 
and 

a = a. + ġa 
Therefore, 

co 1. 

aM AM 
and 

doe l 

da ^a 


The chain rule for partial derivatives yields 


aD _ 3D y 
y (M_+0AM,a_+ódo) = zy Ma) 


Taking the partial derivative of equation (C.19) with 


respect to 0, we obtain 
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dD(M, a) 
98 


20 
3M 


(C.20) 


(0:235) 


(Crer) 


(0.23) 


(C.24) 





aD TUS a Bei] cu 
aw (Ma) = z l Ê DMa) Ê DM pan) + E DM ypa) 


20-1 
+ Ë D(M -$ DO. ou) * S5 D(, 1a.) — (C.25) 


- 20 D(M, sa.) ] 


Using similar procedures, we may derive the remaining 


expressions, 





3D ll 8 20-1 0 
TA (M,a) = BE [ T D(M. 15a. 1) += D(M, ,a.. 4) =i D(M. 4,0. i) 
0 26+1 0 

Pag Elcg pen ze Duis) me aoc (CE 

= 20 DM, 0,3) 
2D (M,a) = EDO ) - 2D(M ) + D(M dl OLD 
aMe 20 (amé s+1* s s* s gs ° 
3^D 1 
T (M,a) = > LD(M, ‚a,_7) = 2D(M , o...) + D(M, 0.472) (C.28) 
da (^a) 
2D y NS l  [D(M ) - D(M ) * D(M ) 
Mag ` SF LAMAa E Sl st1?"s-1 st1 st] 

(C .29) 
- D(M. 193,47) 
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APPENDIX D 


EMPIRICAL RELATIONS 


In the problem treated in Section VI empirical relations 


are used for 


o = f(H), (D.1) 
a = f(H), (D.2) 

and 
M, = £CH) (D63) 


where the parameters are air density ratio (o), speed of 
sound (a), maximum Mach number (My? > and altitude (H). The 


air density ratio is defined as 


ll > 


G 2 (D.4) 
Po 

where o is sea level standard day density. This Appendix 
presents these empirical relations and compares the values 
obtained from these relations with standard atmospheric 


conditions. 


l. Air Density Ratio 


The empirical relation used for air density ratio is 


-c,h 
g = e + cnh (DAS) 


17 





where 


H 
h = Ar (D.6) 
and 
H, = 36,089 ft. (Dep 
e, = 1.54100 (D.8) 
Cy = 1.80445 (D.9) 
C2 = 0.4130 (D.10) 


Figure 50 is a plot of the values of o obtained from 
equation (D.5) compared to those obtained from standard 


atmospheric tables. 


* 


2. Speed of Sound 


The empirical relation used for the speed of sound is 


e 
i 


= a (1-coh) "nh «-I DoD) 


971 ft./sec. , h» 1 (D.12) 


where the parameter a, is the speed of sound at sea level 


ona standard day; that is 


1116.89 ft-/sec. (D.13) 


f 
i 


and 


0.1331 (D.14) 


O 
l! 


Figure 51 is a plot of a vs. H. The expressions (D.11) and 
(D.12) duplicate exactly the values obtained from standard 


atmospheric tables. 
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air density (slugs/ft.? X 1073) 


24 


2.0 


= 
o 


= 
N 
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— — standard day 
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Air density vs. altitude 
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speed of sound (ft. /sec.) 


1100 


1050 

1000 

950 

900 
EXE 0 01 — y L4 
0 20 40 60 80 100 


altitude (ft. x 102) 


Figure 51 


Speed of sound vs. altitude 
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3. Maximum Mach Number 

The empirical relation used for the maximum Mach 
number of the aircraft for the problem solved in Section VI 
is 


-2.lh 


Nw ee Te (D. 15) 


Figure 52 is a plot of equation (E.9) along with the actual 


restrictions of the aircraft under consideration. 


2.2 


2.0 





actual 
—— — empirical formula 





1.6 


14 


Placard Mach number 


1.2 


1.0 
0 20 40 60 3 80 100 
altitude (ft. x 10") 
Figure 52 


Placard Mach number vs. altitude 
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APPENDIX E 


A CONVEXITY THEOREM 


In this Appendix the following theorem on convexity is 
proved. 
Theorem 1: If f(x) is convex on R” where ER”, and f > 0, 
then es is convex on R” where K is any positive integer. 
This theorem is proved by mathematical induction. First, 
the following theorem is proved. 
Theorem 2: If f(x) is convex on R! where a and f » O, 
then f (x) is convex on RÜ, 


Proof: The function f(x) is convex if 


f[Ax, * (1-3)xi] s Af(x54) * (1—))f(xi) (E.1) 
for all X45» X5 and Ae [0,1]. Squaring both sides of 


inequality (E.1) we have 
Felix, + (1-30x,1 « Df (x) * (1-23)£ (4,02? ..— (E.2) 


EncgccuscemonEtboenmedgnality (6.2) is retained as f > 0. To 


prove that f*(x) is convex, it must be shown that 
2 2 2 
DP CLEMENT Gt G-»Xf (x) . (E.3) 


Observing inequality (E.2), (E.3) is seen to be a true 


inequality if it can be shown that 
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[Af(x4) + Q-2)fG421* « AL x,) + (1-1)£ (x,) - (R.N) 


To show that inequality (E.4) is true, we proceed as follows. 


Since) e [0 1], it;is truesthat 
JE ea «ops ) - p(x 512 (E.5) 
22 su. et i 
Expanding inequality (E.5), we have 
A[f(x., - f(x.)]“ < f (x.) - 2f(x.)f(x,) + f2(x.) .(E.6) 
ae a c. Eo SO ea Eu. s. 
Rearranging inequality (E.6), we have 
Mee) meee * 2r(x30f.) - f^(x) « f^(, (E: 7) 
e Ss m E si = mot i 
Subtracting DEKO from both sides, we obtain 


Af(x4) - f(x,)1^ * 2f(x,)[f(x4) - fG23 « f*() - f*03) 
(E.8) 


Multiplying inequality (E.8) by A , we have 


A f*(x4) - 2A7f(x,)f(x,) * A f" (x,)  2Af (x, )£ co) (E.9) 


2 2 2 
= aa (x) < Af (x5) - Af (x) 
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Adding fm, ) to both sides of inequality (E.9), we obtain 


A*f^(x,) - 2A f) fO) € 2Af Gg £9) + f° (a) 


(ELO) 
2 2.2 2 2 2 
-2\f (xi) + AL (x,) SNR (x5) = Af (xi) + f (xi) 
Simplifying inequality (E.10), we obtain 
2 2 2 
[Arlx,) + IE I] «s Af") * OG) . (GM 


This is the inequality we set out to show. The theorem is 


proved. 


Second, the following theorem is proved. 
Theorem 3: If r (x) is convex on R? where x eR”, and f > 0, 


Kee 


then f (x) is convex on Ro. 


Proof: It has already been shown that if f(x) is convex and 


f > 0, then f° (x) is convex. It is now assumed that 
K K K 
f [Xx, + (1-A)x,] < Af (xy) + (1-A)f (xj) . CHI 


Multiplying both sides of inequality (E.12) by f[Ax, + (1-A)x, ] 


a positive quantity, we obtain 


£p € -23)542 « Et Go) Q1) £^ G3 ) 1G Dax Qi) x, 1) 
(E.13) 
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Substituting the expression (E.1) into inequality (E.13), 
we have 
Al BR RIIIE Fa. 
2 2x: epos oe x e AT 
(E.1!) 


To prove that gt ex) is convex it must be shown that 


nor Kl x 


[ax,*(1-3)x4]. « a£ (5) e nr (a) .— (E15) 


Observing inequality (E.14), (E.15) is seen to be a true 


inequality if it can be shown that 


[Xf (x,)#(1-A) £9 (x, )JEAL(x5)+(1=A)£ (24) < at h(x) + ag) 
(E.16) 


To show that inequality (E.16) is true we proceed as follows. 


The expression 
pic =] f(x. IPG.) - f(x.)J (E.17) 
xd ~l ~2 ~l j . 
ls always a positive number because the signs of the expres- 
sions in parentheses in expression (E.17) must be the same. 


The following inequality holds 


ALf (x) - t (x )JEf(x5)-£ 6402. € Ef Go) 5 Gu Ee Go) (x) 1. 
PERIOD 
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Expanding and multiplying by \, we obtain 


NL NL ra (x4 )f(x4)-f(x4 um (x uu a d < 
(B19) 
E EEE CTO 
Rearranging inequality (E.19) and subtracting DE) 
from both sides, we obtain 
A E£P T (x, )- nA (c )f Go) a AAA 
(E.20) 


+A RIA P(e 2 (x1 < ALLE x) 0, 


Further expansion and rearranging yields 


Af ee en (x,)f(x.)-A ep (x5) £ (x5) FAL Ox, Tus (x.)-À "f(x. y AE (x,) 


(E 21) 
ERIT KFI 


a RE eg) SAT TRIAL ) 


(xi 


K+1 


Adding f (xi) to both sides and rearranging further, we 


have 


ph tl 


MEL AMIA E a, (HACIA f Gc) f Go) (i) * e H (xo) 
<A AL )+(1- V (x. ) (B-22) 
or 
[Af (x5)+(1-A) f(x) ITAL (X5)#(1-A) E(x) < APES) # (LAE (x) 
(E.23) 
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This is the inequality we set out to show in (E.16). The 


theorem is proved. 


It has been shown that if f(x) 1s convex and f « 0 then 
f^ (x) is convex. By Theorem 3, £3 (x) ms convex. By 
Theorem 3 again, f(x) is convex. This reasoning can be 
followed for all powers K where K is a positive integer. 


The basic theorem is established. 
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